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Synopsis 


Flow and heat transfer characteristics in a rectangular channel, containing built- 
in vortex generators of both the slender delta-wing and winglet-pair type have been 
analyzed by means of solution of the full Navier-Stokes and energy equations. Each 
wing or winglet pair induces the creation of streamwise longitudinal voitices beliind it. 
The spiraling flow of these vortices serves to entrain fluid from their outside into their 
core. These vortices also disrupt the growth of the thermal boundary layer and serve 
ultimately to bring about the enhancement of heat transfer between the fluid and the 
channel walls. The geometrical configuration considered in the study are representative 
of single elements of either a compact gas-liquid fin-tube crossflow heat exchanger or a 
plate-fin crossflow heat exchanger. Physically, these vortex generators can be mounted 
on the flat surfaces of the above mentioned heat exchangers by punching or embossing 
the flat surfaces. They can also act as spacers for the plate-fins. Because of the favorable 
pressure gradient in the channel, the longitudinal vortices are stable and their influence 
persists over an area which is many times the area of the slender vortex generators. 

The present study gives a quantitative performance data for a delta wing and 
winglet-pair in a channel with different angles of attack and a wide range of operating 
conditions. The study is also extended for turbulent flows where the closure problem is 
solved with the well known k — e model. 

The time averaged Navier Stokes equations for incompressible flows where 
the Reynolds stresses are expressed via the eddy-viscosity concept have been used as 
governing equations. These equations are written in a Cartesian tensor form as 
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where Uj and 9 are nondimensional time-mean velocity components and temperature, 
is the nondimensional kinematic viscosity, is the nondimensional eddy diffusivity 
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and fcji is the nondinaensional kinetic energy of turbulence. The subscripts i and j can 
take the values 1, 2 or 3 in the three coordinate directions Xi, Xo and X3 respectively. 
The Xi, X% and X^ are equivalent to X, y and Z in Cartesian coordinates. The 
Reynolds number i2e is defined on the basis of average axial velocity at the inlet and 
the laminar viscosity as Re = tJoH jv. The Prandtl number of the fluid is Pi . 

The turbulent kinematic viscosity vt,i is given by 

ut,n=C^Rekll€n ( 4 ) 


In the above equation is a constant (equal to 0.09) and e„ is the dissipation rate of 
turbulent kinetic energy (normalized with respect to Uq/H). The nondimensional eddy 
diffusivity at^n is given by 

at,n = C^RePrkl/{at e„) (5) 

where at is the turbulent Prandtl number. Generally, <7^=0. 9. The modeled transport 
equations for and €„ are 
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Here G„ is the generation of turbulent kinetic energy(nondimensioaal) given by 


dUX dU. 

Re [dXj dXi) dXj ^ ’ 

and crfc Cl, C2 are empirical constants : ak = 1.0, erg = 1.3, Ci = 1.44, C2 = 1.92. 
These are based on wide range of experimental data. 

The equations are solved by a numerical procedure based on MAC (Marker and 
Cell) method. It may be mentioned that for the laminar flow computations, ut^n and A;,^ 
in equation (2) and at^n in equation (3) are taken as zero. Possibly it is needless to say 
that for the laminar flows, the barred quantities should' be read as unb<irred variables. 

Figure S-1 compares the numerical flow visualization with its experimental counter- 
part due to a research group in the Ruhr University Bochum of Germany. The as[)ect 
ratio, A, and the angle of attack, l3, of the winglet for this study are 1 and 35'^ respec- 
tively. In the experimental study, at the rear edge of the winglet, a primary longitudinal 
vortex of elliptical shape, is seen. The Reynolds number of this flow visualization study 
is 1350 and the flow is laminar. In the nmnerical flow visualization given, the secondary 
velocity vectors enable a view of the vortex traiispcn-t at different regions (d’ the cross- 
sections of the channel. A qualitative comparison between tlu^se numerical and the 
experimental flow visualizations implicates favorable agreement b(!twtien the two. How- 






ever, a detailed evaluatioa of the performance parameters with regard to the enhance- 
ment of heat transfer using various wing-type vortex generators has been accomplished 
in this study. 

From heat transfer point of view, delta wing is found to be more effective than 
winglet-pair. However, in the case of convective heat transfer two type of losses of useful 
work is encountered; one due to fluid friction resulting in viscous heat generation and 
other due to heat transfer across finite temperature gradient. Both these phenomena 
are manifestations of irreversibility and therefore for a complete evaluation of a heat 
transfer augmentation technique a detailed thermodynamic analysis is required. Thus, 
the present study incorporates an in-depth discussion about the influence of vortex 
generators (wings /winglets) on entropy generation and irreversibility. 

It may be mentioned that the study of laminar flow in the present work is not for 
computational simplification. Usually, the fin spacing is so small and the mean velocity 
range is such that the flows are often laminar. However, for some special applications 
the velocity becomes very high and one may encounter turbulent flows in plate-fin heat 
exchangers. In the present study, results for turbulent flow have been obtained for a 
range of Reynolds number between 5000 and 100000. An effort has been undertaken to 
evaluate the accuracy of A: — e model in predicting such flows. In order to accomplish 
this, our computational results have been compared with a well documented experiment 
of another research group in Stanford University of the USA. The computational results 
agree reasonably well with their experimental counterpart. 
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Chapter 1 


Introduction 


1.1 Description of the Problem 

Fill-tubes are generally used in gas-liquid heat exchangers where the gas is in the 
crossflow to the tubes and the liquid flows inside the tubes. Figure 1(a) shows the 
arrangement of the core region of a fin-tube crossflow heat exchanger. The purpose of 
the fin is to enhance the heat transfer coefficient on the gas side to a value comparable 
to the liquid side. The area ratio of the fins to the tubes varies depending on the 
application and may reach a value of 30. Plate-fin heat exchangers also find application 
in many industrial processes. Such heat exchangers may be used to exchange energy 
between two different gases or liquids. Figure 1(b) illustrates the basic arrangement 
of the core region of a plate-fin heat exchanger. In forced convection heat transfer 
between a gas and a liquid, the heat transfer coefficient of the gas may be 10 to 50 
times smaller than that of the liquid. In order to increase the heat transfer on the gas 
side, the fin geometry can be manipulated (Webb, 1987) in various ways. In addition, 
compactness of the heat exchangers is demanded for effective utilization of energy and to 
bring about a save on material and space. The performance of heat exchanger surfaces 
can be enhanced by moimting protrusions on the surfaces. The offset strip-fin and the 
louvered fins have been widely used for this purpose (Webb, 1987). The basic mechanism 
involved in the above mentioned cases is the periodic interruption of the boundary layer 
by separation, wake recovery and generation of thin developing boundary layer with 
high Nusselt numbers. 

A somewhat different concept for the reduction of the thermal resistance is to induce 
longitudinal streamwise vortices in the flow field. The longitudinal vortices can be 
generated by mounting delta wing or winglet-type vortex generators on the flat surfaces 
(Fig. 1.2). The longitudinal vortices develop along the side edge of the wing-shaped 
vortex generator due to the pressure difference between the fi'ont surface facing flow and 
the back surface. These streamwise vortices interact with an otherwise two-dimensional 
boundary layer and produce a three dimensional swirling flow that mixes near wall 
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Liquid 



(a) 


Figure 1.1: Typical arrangement of heat exchanger cores (a) Gas-liquid fin-tube 
cross-flow (b) Plate-fin (single or multi pass) 

fluid with the free stream. This mechanism strongly enhances the exchange of fluid 
between the wall and the core region of the flow field which causes high heat transfer 
augmentation. The additional pressure losses are modest because the form drag for 
such wing-t 3 rpe slender bodies is low. A complete numerical evaluation of heat transfer 
in a fin-tube or plate-fin crossflow heat exchanger would be extremely difficult. In 
order to understand the basic mechanisms involved in enhancement of heat transfer 
and to evaluate the performance parameters quantitatively, a need is felt to analyze 
three-dimensional flow and heat transfer in a horizontal channel with built-in vortex 
generators. 

A number of experimental investigations (Edwards and Alker, 1974; Russel et o/., 
1982; Fiebig et ai, 1986) in fields relevant to the present problem have been reported 
in the literature. Computational studies on related topics have been performed by 
Fiebig et al. (1989) and Biswas and Chattopadhyay (1992) for laminar flows. Arising 
out of the practical and academic importance of the configmation of interest, it is 
conjectured that a detailed analysis of the longitudinal vortices embedded in a shear 
flow will provide a better understanding of the complex phenomena and facilitate a 
more accurate estimation of the pertinent performance parameters. 


1.2 The Scope of the Present Work 

Since, the details of the interactic^us between the longitudinal strt'arnwise vortices, the 
enhancement of heat transfer and the associated flow los.s<^s for laminar and turbulcuit 
flows are not still clearly understood, there is a general agre('meut that more rest'arch 
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Figure 1.2: Protrusions in the form of (a) delta wing and (b) winglet pair on the flat 
surface to enhance heat transfer 

on the subject is required. 

In this thesis, a detailed three-dimensional numerical model has been formulated 
which is able to provide a clear understanding of the physics of the flow. Grid indepen- 
dent numerical solutions have been envisaged which carefully track the development of 
swirling flows, in order to estimate more accurately the performance parameters. In the 
numerical model, the full N avier-St oke s equations together with the governing equation 
for energy have been solved in a rectangular channel with built-in wing-type vortex 
generators. The present investigation also highlights the thermodynamic optimization 
of the heat transfer with respect to variation in geometrical parameters involved in 
transport-enhancement. An attempt has also been made to compute turbulent flows 
with suitably adapted k — e model. An assessment regarding the validation of the model 
has been tried as well through the comparison with experimental results due to Pauley 
and Eaton (1988a, b). 


1.3 Layout of the Thesis 

Chapter-2 of the thesis provides a review of literature relevant to the basic mech- 
anisms involved in augmentation of heat transfer with special emphasis on generation 
of longitudinal streamwise vortices. A concise review of the modeling aspects of tur- 
bulent flow has also been focussed in this chapter. The mathematical formulation of 
the problem for simulation is presented in chapter-3. In this chapter, the geometry, 
governing equations, boundary conditions, method of solution of the governing equa- 
tions to obtain the nmnerical results are described in detail. Chapter-4 discusses the 
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results for laminar flows. Basically a detailed appraisal of the performance parameters 
is accomplished in this chapter. Two different wing-type vortex generators have been 
investigated numerically, viz., delta wings and delta winglet pairs. In chapter-3, relative 
performance of the wing and winglet-pair has been compared. The basis of compari- 
son has been a second law analysis. Most convective heat transfer processes encounter 
two types of losses, namely, losses due to fluid friction and those due to heat transfer 
across a finite temperature gradient. These two phenomena are the manifestations of 
thermodynamic irreversibility and the investigation of a process from this standpoint is 
known as second law analysis. Chapter-6 presents the results for turbulent flows. Tliis 
chapter serves bivalent purposes. On one hand, a critical appraisal of the numerical 
model is accomplished through comparison with a well documented experimental data 
base of a research group in Stanford University and on the other hand, a parametric 
study is presented predicting the effects of various geometric and operating variables 
in the turbulent regime. Chapter-7 includes the concluding remarks and the scope for 
further research. 



Chapter 2 

Review of Literature 


2.1 Introduction 

The objectives of the present study have beeti briefly outlined in the previous chapter. 
A variety of experimental and a considerable amount of analytical and computational 
research have been carried out on the enhancement of heat transfer. Especially, the 
enhancement of heat transfer through manipulation of surface geometry has concerned 
many researchers and practitioners since the earliest documented studies of heat trans- 
fer. In this chapter a survey of the relevant literature is presented to indicate the extent 
of work already reported in open literature pertaining to the enhancement of heat trans- 
fer by introducing protrusions mounted on the heat transfer surfaces. In order to analyze 
the flow structure and heat transfer in such applications a detailed computational study 
is needed. The flow regime for such applications can be both laminar and turbulent. In 
the present Uterature review the attention is also focussed, by and large, on important 
numerical techniques developed in the recent past for solving the governing conserva- 
tion equations for complex flows and heat transfer. Although the immerical solution 
of the governing conservation equations for laminar flows has become a realizable goal, 
the solution for turbulent flows even for simpler geometry is indeed a formidable task. 
The relevant literature on turbulence modeling has also been looked into with intent to 
understand various modeUng and computational strategies. However, this stu'vey also 
helps in suggesting further work that should be carried out to accomplish the objectives 
entunerated in the earlier chapter. The literature is reviewed from three different view 
points. In the first section, an overview of past work done by various researchers in the 
area of augmentation of heat transfer is presented and discussed. The second section 
presents a review of investigations pertaining to the numerical methods with regard to 
the task of computing flow fields in complex geometries. Basically, this section deals 
with the available schemes for solving the full Navier-Stokes equations. In the third, a 
comprehensive review on modeling aspects of turbulent flow has been midertaken with 
a view to suggest a suitable turbulent-closure-model for our problem. 
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2,2 Enhancement of Heat Transfer by Mounting 
Protrusions on the Surfaces 

The subject of enhanced heat transfer is of serious interest in heat exchanger appli- 
cations. Concepts like space and energy optimization in power and process industries 
call for more and more compact designs of heat exchangers. Specially designed surfaces 
use protrusions mounted on them to enhance heat transfer. Some typical examples 
of the enhanced heat transfer surfaces which are very popular in diflFerent industrial 
applications (Webb, 1987) are shown in Fig. 2.1. There have been a number of re- 
cent survey articles and handbook sections prepared that deal with the augmentation 
of heat transfer for different applications (Bergles, 1978, 1983, 1985). Special surface 
geometries bring about the transport-enhancement by establistiing a higher "hA" per 
unit base surface area. It may be mentioned that (1/hA) is the thermal resistance as- 
sociated with heat transfer by convection at a surface. Some well known methods of 
manipulation of surface-geometry are; 


(a) Inserting twisted tapes or turbulence promoters for internal flows in circular 
tubes. 

(b) Roughening the surfaces for channel flows in closely spaced parallel plate 
channels, typical of plate -fin and plate type heat exchangers. 

(c) Improving flow structure in external flows normal to tubes and tube banks. 


2.2.1 Internal Flows in Tubes and Ducts 

Some useful methods of enhancing heat transfer on the inner surfaces of the tubes 
have been documented by Webb (1987). For tubular heat exchangers, enhancement is 
desirable for the inner side, outer side or both sides of the tube. Fig. 2.2 shows some 
doubly enhanced tubes where the enhancement is applied to both iumu' and outer sur- 
faces. Generally the enhancement applied to the innerside is caused by imparting a swirl 
to the flowing fluid. Date (1974) h<xs performed numerical predictions for the twisted 
tape in fully developed laminar flow for constant heat flux with constant properties. His 
analysis shows that the Nusselt number increases with increasing Prandtl and Reyiudds 
numbers, contrary to laminar flow in smooth tubes which is independent of Re and 

Pr. Date (1973) has also proposed a correlation for the friction factor for the range 
5000 < Re < 70,000. 
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Figure 2.1; (a) Plate-fin heat exchanger and its surface geometries : (b) plain rectangular 
fins (c) plain triangular fins (d) wavy fins (e) offset strip fins (f) perforated fins (g) 
louvered fins; after Webb (1987) 


Bergles et al. (1980) have reported a bibliography of different augmentation tech- 
niques. Nondecaying swirl flow has been the subject of many investigations (see Kreith 
and Margdis, 1959; Gambill and Bundy, 1963; Sniithberg and Landis, 1964; Hong and 
Bergles, 1976). Blum and Ohver (1966) and Migay and Golubev (1970) have shown 
that there is significant increase in heat transfer due to free swirling flow. An analytical 
study of the heat transfer characteristics in decaying turbulent swirl flow generated by 
short twisted tapes placed at the entrance of the test section is carried out by Algifri and 
Bharadwaj (1985). They have shown that the augmentation in the local heat transfer 
can be as high as 80 percent and the augmentation is worthwhile up to an initial length 
of about 60 tube diameters. Bergles and Joshi (1983) provide an excellent summary 
of performance data on the twisted tape-inserts for the constant wall temperature and 
constant wall heat flux boundary conditions. Sparrow and Chaboki (1984) have per- 
formed an experimental study on swirl affected turbulent air flow and heat transfer in 
a tube. The swirling motion has been found to enhance heat transfer substantially in 
the initial portion of the tube. Compared with the heat transfer coefficient encountered 
in the conventional thermal entrance region in a nonswirling pipe flow, those associated 
with swirl are found to be remarkably greater. Junkhan et al. (1985) have conducted 
experimental studies of three different type of tape inserts for fire tube boilers. These 
tapes are termed as turbulators. Two commercial turbulators consisting of narrow, thin 
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(b) 


Twist ratio 
= rt = p/adi 



= Axial distance through 
which the tape is twisted 
by 180 “ 


Figure 2.2: Various metliods of making eiihcxnced. tubes: (a) iielical rib on inner surface 
and integral fin on outer surface (b) internal fins on inner surface and porous coating on 
outer surface (c) twisted tape insert on inner surface and integral fins on out(U' surtacc' 
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metal strips bent and twisted in zig-zag fasliion to allow a periodic contact with tube 
wall, have displayed 135 and 175 percent increase in heat transfer coefficient at a liigh 
Reynolds number. A tliird turbulator consisting of a twisted strip, with width slightly 
less than the tube diameter, has provided a 65 percent increase in the heat transfer 
coefficient while the increase in the friction factor is small as compared to the other two 
turbulator inserts. 


2.2.2 Flows in P 2 irallel Plate Channels Formed by Two Neigh- 
boring Fins of Plate -fin or Fin-tube Heat Exchangers 

It may be mentioned here that we are particularly interested in enhancement of heat 
transfer in the gas side of fin-tube heat exchangers (with fiat fins) and enhancement 
of heat transfer in plate-fin heat exchangers. With this intent, we would like to study 
different investigations related to augmentation of heat transfer suitable for preceding 
applications. 

Augmentation of heat transfer is of special interest in channel flows where the rate of 
the heat transfer between the fluid and the channel walls deteriorates as the boundary 
layer grows on the channel walls and the flow tends to become fully developed. Pro- 
trusions can be mounted on these channel walls in order to disrupt the growth of the 
boundary layer and thereby enhance the heat transfer between the flowing fluid and the 
channel walls. Two relevant applications using this kind of flow configuration are the 
heat transfer between the gas and the fin in the case of gas-liquid fin-tube cross-flow 
heat exchangers and the heat transfer between flowing fluid and plates in the case of 
plate-fin heat exchangers (Fig. 1.1). The evolution towards a fully developed flow can 
be disturbed by using a multi-louvered surface geometry for the plates. Investigation 
by Achaichia and Cowell (1988) provides a detailed performance data for louvered fin 
surfaces. However in using louvered fins, enhancement is obtained at the price of high 
pressure drop. To circumvent this difficulty, protrusions in the form of slender delta 
wings or winglets can be deployed (Fig. 1.2). As shown, the base of the wing remains 
attached to the fin and the apex faces the incoming stream with an angle of attack with 
this configuration. The longitudinal vortices are generated along the side edge of the 
wing-shaped vortex generator due to the pressure difference between the front surface 
facing the flow and the back surface (Fig. 2.3). These longitudinal vortices, generated 
by the vortex generators, can be made to disrupt the growth of boundary layer in a 
channel by exchanging the fluid fi:om the near-wall-region with the channel-core-region 
and thus they can serve to enhance the heat transfer rate while producing less of a 
pressure drop. Use of longitudinal vortices for boundary layer control is well known 
(Pearcey, 1961) and the vortex generators are used in commercial air-planes for this 
purpose. 

Formation of streamwise longitudinal vortices behind a slender aerodynamic object 
is a research topic of considerable interest for several years. Both the theoretical and 
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Figure 2.3: Vortices formed by a slender delta wing; after van Dyke (1982) 


experimental investigations on flow past a delta wing have been conducted and reported 
in literature by a number of researchers. Among them Winter (1956), Fink (1956), 
Marsden et al. (1958) and Lawford (1964) are noteworthy. The physics of vortex 
formation has been furnished by Brown et al. (1954), Mangier et al. (1959) and Smith 
et al. (1968). Use of thin aerofoil theory provides good approximation of the strength of 
the leading edge vortices and the lift, though in trailing edge part their results show some 
discrepancy with experiments. Rehbach (1976) and Kandil et al. (1976) have obtained 
better results with non-conical theory. Computational solution of non-conical flow field 
and pressure distribution on wing surface has also been attempted by Johnson et al. 
(1976) and Weber et al. (1976). Their study suffers from some numerical complications 
in the vicinity of the trailing edge. Hummel and Sriuivasan (1967) have made important 
contribution in revealing the complex flow structure behind a delta wing. They hcive 
conducted experiments and presented pressure distribution and vortex structure of the 
flow around a delta wing of unit aspect ratio with an angle of attack of 20". 

Thomas et al. (1990) have computed low-speed laminar flow over a low aspect ratio 
delta wing up to 40'' angle of attack using an upwind biased flnite volume algoritlun. 
The differ enciiig scheme used are second-order accurate and a inultigrid algorithm is 
('mployed to promote couvergenc<i to steady state. The predicted results and lift c:o('ffi- 
cient of the wing have remarkable agreement with expr'riments due, to Hummel (1973). 
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At 40" angles of attack, a bubble-type vortex breakdown is ('vident in the computations. 

Ex]XTimental investigations due to Edwards and Alker (1974), Russels et al. (1982), 
Fiebig et al (1986), Fiebig et al (1991) and Tiggelbeck et al. (1992) can be reh^rred 
to in connection with augmentation of heat transfer by means of longitudinal vortices. 
Experimental investigation due to Fiebig et al. (1991) is the first systematic study to 
compare the performance of different kinds of vortex generators, viz., delta wing, rect- 
angular wing, delta winglet pair and rectangular winglet pair in the Reynolds number 
range of 1360 and 2270. Their observations depict that the delta wing is the best vortex 
generator from heat transfer point of view. Another important feature of their obser- 
vations is that the heat transfer coefficient increases with increase in angle of attack till 
the vortex break down takes place. Tiggelbeck et al. (1992) use multiple rows of vortex 
generators in an aligned arrangement witliin a channel and observe their influence on 
flow structure. The flow structure in the wake of the second row is qualitatively similar 
to that of the first row. Their flow visualization by a laser light sheet technique has 
revealed that the concentrated vortex pair generated by a small aspect ratio delta wing 
at large angle of attack has elliptic shape due to the influence of channel walls. They 
also observe that the peak value of spanwise-averaged Nusselt number at the wake of 
the second row is strongly dependent on the spacing of the two rows. 

Computational studies on related topics have been performed by Fiebig et al. (1989) 
and Biswas and Chattopadhyay (1992) for laminar flows in a geometrical configuration 
of delta wing placed inside a channel. Both studies have discussed the influence of angle 
of attack and Reynolds number on velocity and temperature fields. 

When the flow velocity is not very high and the temperature difference between 
the body surface and ambient fluid is large, flow and heat transfer characteristics are 
strongly influenced by thermal buoyancy forces. Even without a vortex generator, buoy- 
ancy can induce longitudinal vortices and multiple plumes which can enhance heat trans- 
fer and alter the length of entrance region in mixed convection flows in ducts (Patankar 
et al., 1978; Incropera and Schutt, 1985; Incropera et al. , 1987). Biswas et al. (1989) 
characterize the mixed-convection condition by buoyancy driven secondary flows that 
form counter rotating longitudinal vortices and enhance heat transfer. In another study, 
Biswas et al. (1990) observe that the critical Reynolds number in a channel flow with 
built-in obstacle may be lower in the case of mixed convection than in the case of forced 
convection alone. It is observed that for a given Reynolds number, the heat transfer 
is improved up to a certain Grashof nmnber and deteriorates if the Grashof number is 
further increased. Nandakumar et al. (1985) have obtained dual solutions for the case 
of a duct flow. In one solution a two vortex pattern emerges while in the other solution 
a four vortex pattern evolves. 

It has already been discussed that the gas side heat transfer coefficient in gas-hquid 
fin-tube crossflow heat exchangex's is small as compared to that of the liquid side. The 
mechanism of heat transfer between the gas and the solid surfaces in such cases is to be 
understood in detail. Let us consider only one tube on a plate-fin (Fig. 2.4). The flow 
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Plate 


Tube 



vortex system 


Figure 2.4: Flow around a tube on a plate 


field consists of a horseshoe vortex system, a dead water zone at the juncture of the tube 
and the plate and a von-Karman vortex street in the middle. In order to enhance heat 
transfer in such a flow configuration, vortex generators can be mounted on the plates 
which makes the flow field extremely complicated. Recently, Dong (1989) and Valencia 
(1992) have conducted experimental investigations to observe the influence of winglet- 
type vortex generators in a channel with a built-in circular tube. Sanchez et al. (1989) 
have computed laminar flows around a circular cylinder in a rectangular channel with 
a pair of delta winglets on the bottom plate of the channel. Their results show that the 
longitudinal vortices generated by the winglets placed in the wake, control the spread 
of the wake zone behind the cylinder and damp the periodic vortex street. Biswas et 
al. (1994) have observed that in the absence of the winglet-type vortex generators, 
relatively little heat transfer takes place in the downstream of the circular tube which is 
a recirculation region with low velocity fluid. However, they observe an enhancement of 
heat transfer as high as 240 percent in the wake region behind the cylinder in presence 
of winglet-type longitudinal vortex generators. 

Recently, computational studies have been extended to periodic configurations wliich 
are realistic in terms of plate- fin heat exchangers with perforated fins or slotted chan- 
nels. Heat transfer surfaces periodically interrupted along the streamwise direction have 
een widely used in order to obtain improved performance of the heat-exchange devices. 
Such an arrangement may be viewed cis a succession of collinear plate-segments aligned 
parallel to the flow with slots between successive plates. Each slot enables the inter- 
ruption of the thermal boundary hiyer where the liiglu'st resistance to the heat flux 
occiu-s and a new boundary layer with lower thermal r^jsistauce begins. It also allows 
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a v('.ry good mixing of tlio flow with .sclf-s’ustaiiK'd oscillation when the' system is oi)er- 
ating above tlu' critical conditions (Reynolds number). The channels wiiich have such 
operating conditions liave been cla.ssified as communicating channels. As it has been 
told, above a critical value of the Reym)lds number (much within the laminar regime), 
the initially steady flow becomes unstable to small disturbances. At this instant, the 
flow is set to have undergone a bifurcation from steady state to a time periodic, self- 
sustaining oscillatory regime. These flows are characterized by supercritical bifurcation 
and narrow band excitation (Ghaddar et al. , 1986a; Ghaddar et al. , 1986b; Amon 
and Patera, 1989). Appreciable transport augmentation takes place in this range of 
Reynolds number (Majumdar and Amon, 1992). 

Webb and Ramadhyani (1985) have analyzed fluid flow and heat transfer characteris- 
tics through a parallel plate channel with transverse ribs which act as vortex generators. 
Significant heat transfer augmentation is obtained in their study for fluids such as water 
or flourocarbons. Acharya et al. (1991) explain vortex interaction between two such 
ribs in a channel for the subharmonic excitation of shear layer. In a similar study, sig- 
nificant augmentation is obtained by positioning vortex generators at key locations in 
the flow ( Myrum et al., 1992). 

Many researchers have observed longitudinal vortices in various complex flow config- 
urations. Naturally, the interaction of such vortices with boundary layer and its effect 
on heat transfer is a subject of interest to them. Some examples are Taylor-Gortler 
vortices in boundary layers on concave curved surfaces, horseshoe vortices formed by 
an obstruction protruding from a surface and wingtip vortices impinging on a down- 
stream surface. The embedded vortex is capable of strongly perturbing the boundary 
layer thickness and influencing the heat transfer characteristics. In addition, longitu- 
dinal vortices usually maintain their coherence over a long streamwise distance. As a 
consequence, the heat transfer effects behind a vortex generator are very persistent. 
Eibeck and Eaton (1987), Westphal and Mehta (1987) have worked extensively in this 
field with more focus to the turbulent boundciry layer. Eibeck and Eaton (1987) have 
conducted experiments on longitudinal vortices embedded in a turbulent boundary layer 
and resultant heat transfer effects. Longitudinal vortices are found to influence the heat 
transfer behaviour significantly. Local Stanton number is increased by as much as 24 
percent resulting in a net increase in spanwise average heat transfer coefficient. Despite 
the presence of turbulent diffusion, the influence of longitudinal vortices on momentum 
and energy transport can be traced at a location as far downstream as 60 wing chords 
behind the delta wing. 

Westphal and Mehta (1987) have studied the effect of oscillating vortex on a turbu- 
lent boundary layer. The meander is simulated by forcing a periodic lateral translation 
of a half-delta wing vortex generator at a low frequency of 1 Hz and one-half amplitude 
of 0.5 cm. The effect of the meander is found to flatten the vorticity contours at the 
stations where they are originally round. The Reynolds stresses are also found to be 
significantly affected. Ligrani et al. (1991) have presented results that illustrate the 
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effects of embedded longitudinal vortices on heat transfer and also the effect of injecting 
a fluid downstream of a row of film-cooling holes in a turbulent boundary layer. In 
some circumstances, the protection provided by film cooling is augmented because of 
secondary flows, which cause extra injectant to accumulate near upwash regions. 


2.3 Development of Numerical Methods for Solving 
Full Navier- Stokes Equations 

2.3.1 Methods Developed on Regular Geometries 

Different solution techniques for full Navier-Stokes equations have been developed 
during the past three decades. As it has been seen, the major difficulty encountered 
during the solution of incompressible flows is the non-availability of any obvious equa- 
tion for pressure. This difficulty can be resolved in the stream function-vorticity ap- 
proach. But the stream function-vorticity approach loses its attractiveness when three- 
dimensional flow is computed because of the absence of a single scaler stream function in 
the three-dimensional space. A three dimensional problem demands a primitive-variable 
approach. Efforts have been made so that two-dimensional as well as three-dimensional 
problems could be computed following a primitive variable approach without encounter- 
ing non-physical wiggles in pressure distribution. As a remedy, it has been suggested to 
employ a different grid for each-dependent variable. Harlow and Welch (1965) have used 
such a staggered grid for the dependent variables in their well known MAC (Marker and 
Cell) method. The MAC method of Harlow and Welch is one of the earliest and widely 
used explicit methods for solving full Navier-Stokes equations. In this method, solutions 
of velocities are obtained in two folds. In the first fold, provisional value of velocity com- 
ponents are computed explicitly using advection, diffusion and pressure gradients of the 
earlier time step. This explicitly advanced provisional velocity field may not ensure a 
divergence free velocity field. In the subsequent second fold, pressure and velocity com- 
ponents are corrected through the solution of a Poisson equation for pressure. A related 
technique developed by Chorin (1967) involves a simultaneous iteration on pressiues 
and velocity components. Vicelli (1971) has shown that the two methods as applied to 
the MAC method are equivalent. The original version of the MAC method luis been 
modified by Harlow and Amsden (1970), Nichols and Hirt (1971) and Hirt and Cook 
(1972) for application to free surface flows. The MAC method uses a layer of irn<iginary 
cells around the boundary of the physical domain necessitate updating of bound<u:y con- 
ditions after every change in interiml velocity and pressure values. The MAC method 
has been extensively used by many researchers to solve flow in complex geometry. For 
example, Braza, Chassaiiig and Ha-Minh (1986) have (obtained the unsteady wake bc'- 
liiud circular cylinder. Mukhoi)adhyay, Biswas and Sundararajan (1993) have obtained 
the periodic wake behind a rectangular obstacle. In fact, the MAC mc^thod has bc'cn 
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successfully used even to simulate liighly unsteady turbulent flows (Robichaux, Tafti 
and Vankfi, 1992). It has bet'ii exi)erienc(!d that the MAC nu’thod is ind(>ed very effi- 
cient in the studies of temporal flow dev(doi)ment but it has stability restrictions on the 
time step values wliich slow down the calculations for steady flow considerably. Since 
implicit methods have no such restrictions, they are obviously very attractive. Patankar 
and Spalding (1972) have introduced an efficient implicit method known as SIMPLE 
(Semi-Implicit Method for Pressure Linked Equations). This method is based on a fi- 
nite volume discretization of the governing equations on a staggered grid. In order to 
improve the convergence involved in the pressure-velocity coupling, several variants of 
SIMPLE algorithm have been developed. The SIMPLER algorithm of Patankar (1981) 
and the SIMPLEC algoritlim of van Doormaal and Raithby (1984) may be referred 
to with regard to such improvement on SIMPLE algorithm. Althougli, the necessary 
changes to incorporate SIMPLEC into a SIMPLE algorithm is minor, the consequences 
can be great as it eliminates the approximations made in SIMPLE while deriving the 
pressure-velocity corrections. A comi:)arative illustration between the operator splitting 
algorithm, viz., the PISO of Issa (1986) and the SIMPLE family of algorithm has been 
reported by Jang, Jetli and Acharya (1986). 


2.3.2 Methods on Complex Geometries 

In regular geometry, the orthogonality of the coordinate frame makes it possible 
to decouple the integrands so as to carry out individual line integration products for 
evaluating the surface integrals. Thus, many researchers have tried to generate orthog- 
onal grids for complex solution domains while considerable effort has also been directed 
towards solving the equations on general non-orthogonal curvilinear grids. 

In order to extend the aforesaid established flow solvers to the complex geometries, 
generalized orthogonal body-conforming coordinate systems have been used in different 
combinations. Pope (1978) has applied the conservative forms of the equations with 
flnite volume method to compute turbulent recirculating flows in diffusers. Gosman 
and Rapley (1980) have used a similar procedure for fully developed flows in ducts of 
arbitrary cross-section. Lawal and Mujumdar (1985) have performed calculations in 
the entrance region of an arbitrary shaped cross section. In this case, collocated carte- 
sian velocity components on orthogonal grids in each cross-sectional plane have been 
solved along with a Poisson equation for pressure. Garg and Maji (1987) have applied 
SIMPLEC method for the solution of viscous flows through periodically converging- 
diverging tubes. In another numerical investigation, Velusamy and Garg (1993) have 
solved the complete set of Navier- Stokes equations for the three-dimensional developing 
flow in elliptical cross-section ducts. 

In non-orthogonal systems, general curvilinear coordinates are employed which ex- 
actly coincide with the boundaries of the physical domain. The governing partial differ- 
ential equations are transformed in terms of these curvilinear coordinates and finite dif- 



16 


CHAPTER 2. REVIEW OF LITERATURE 


ference representations are derived in the transformed domain. Gal-chen and Somerville 
(1975) and Faghri et al. (1984) have used algebraic coordinate transformation. While 
the former use a MAC-like grid arrangement coinciding with the grid lines, latter used 
covariant components of velocities. Thompson et al. (1982) have made remarkable 
contributions to the development of numerical grid generation techniques for solving el- 
liptic partial differential equations related to both external and internal flow problems. 
Vanka et al. (1980), Maliska and Raithby (1984) and Shyy et al. (1985) have extended 
the available solvers of different flow problems to irregular geometries. The basic flow 
solvers have been based on the SIMPLE algorithm. The adaptive grid procedure of 
Acharya and Patankar (1985) for parabolic flows may be mentioned here for its detailed 
discussion on the choice of different discretization schemes and their performance in 
the adaptive grid layout. However, application of finite volume methods using non- 
orthogonal coordinates and collocated grids are reported by Rhie and Chow (1983) and 
Peric (1985). Peric et al. (1988) have shown that the collocated arrangement converges 
faster than the staggered variable arrangement and has advantages when extensions such 
as multigrid techniques and non-orthogonal grids are considered. It may be worthwhile 
to note from Benard and Thompson (1984) that while staggered grid is a direct solution 
to avoid pressure split with conventional finite difference interpolations of flow variables, 
proper interpolation of velocities still remains crucial to avoid oscillatory velocity field 
even on a staggered grid. An appropriate pressure interpolation avoids pressure splits 
in these methods of collocated grids. A finite volume based procediue has been very 
successfully developed by Majumdar et al. (1992) using such collocated velocities and 
pressure. The concept of "momentum interpolation" of cell face pressiures from nodal 
values has been established showing equivalence in the methods of staggered and collo- 
cated arrangements. Mukhopadhyay et al. (1993) have developed a numerical method 
for predicting viscous flows in complex geometries. Integral mciss and momentum con- 
servation equations are deployed and these are discretized into algebraic form through 
niunerical quadrature. The physical domain is divided into a number of non-orthogonal 
control volumes which are isoparametrically mapped on to standard rectangular cells. 
Numerical integration for unsteady momentum equations is performed over such non- 
orthogonal cells. The algorithm is tested on some complex problems. The results exhibit 
good accuracy and justify the applicability of the algorithm. 


2.3.3 Modelling of Convection-Diffusion Equations 

During the development of flow solvers, emphasis is given to efficient handling of non- 
linear convection terms in the flow equations. As such entire treatment of convection- 
diffusion equations is somewhat tricky due to the presence of the first-order nonlinear 
spatial derivatives that describe the convective flux. Especially, in the case of tran- 
sient equations, the interaction between space and time dependcuice in tlui convection- 
diffusion equation requires that the temporal and spatial derivatives be accurately ap- 
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proximated in order to obtain a satisfactory solution. Solutions based on conventional 
finite difference or finite element spatial discretizations suffer from spurious oscillations 
due to space-centered discretization of the convective terms. The possible improvement 
is to use a first-order upwinding technique. However, the first-order upwind schemes 
introduce significant amount of artificial (numerical) diffusion into the physical prob- 
lem. Roache (1972), Roberts and Wiess (1966) and Bozeman and Dalton (1973) have 
criticized the idea of introducing this "artificial viscosity" into the physical problem 
through upwind schemes which are first-order accurate. Roache (1972) has estimated 
the amount of false diffusion due to such schemes and proved that its contribution may 
be too high to be ignored. Spalding (1972) and Cristie(1976) have presented a local 
analytical solution of convection-diffusion equations which is superior to central differ- 
encing schemes, when the local Peclet number of the grid is large. Leonard (1979) has 
developed a second-order upwinding scheme, known as Quadratic Upstream Interpo- 
lation for Convective Kinematics (QUICK) in order to avoid the stability problems of 
central differencing. However, this second-order upwinding involves a dispersive third- 
derivative in the truncation error. Raithby (1976) has introduced two schemes, namely. 
Skew Upstream Differencing Scheme (SUDS) and Skew Upstream Weighted Difference 
Scheme (SUWDS) to reduce the error, brought about in the flow region when the grid 
line and velocity directions are not aligned. These schemes are capable of reducing 
the disadvantages of upstream differencing. A simple stable mass-flow-weighted skew 
upwind scheme has been proposed by Hassan et al. (1983). Leonard (1984) has shown 
that the third-order upwinding provides a robust discretization procedure in the mod- 
eling of convection-diffusion terms of Navier Stokes equations. Runchal (1987) has 
modified central difference scheme by introducing a controlled amount of numerical 
diffusion based on local gradients. This proposed technique for Controlled Numerical 
Diffusion with Internal Feedback (CONDIF) eliminates the undesirable feature of the 
spaced-centered schemes. Shirayama (1992) has constructed the modified third-order 
upwinding schemes for stretched meshes using Lagrangian interpolation formula. This 
method prevents numerical instability, preserves third-order accuracy and produces re- 
alistic results. Although upwind differencing scheme has its dissipative nature, it is 
probably the best choice for high Reynolds number (Re > 1000) flows. Patankar (1980) 
has mentioned about various approximative schemes, namely, upwind, hybrid, exponen- 
tial and power-law schemes and discussed on their relative merits and demerits. Galpin 
and Raithby (1986) have linearized the nonlinear convective terms by Newton-Raphson 
linearization method and found that such a linearization enhances convergence, espe- 
cially when grid curvatiue effects are important. Pollard and Siu (1982) and Thakur 
and Shyy (1993) have successfully used upwind, hybrid and QUICK schemes in compu- 
tation of complex flows. Tamamidis et al. (1993) have made use of hybrid, third-order, 
QUICK and fifth-order upwind schemes for modeling turbulent flows using SIMPLEC 
algorithm. The results show that accurate solutions can be efficiently obtained on grids 
of moderate size by using high-order-accuracy schemes. 
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It is clear that substantial progress has been made on development of higher or- 
der schemes which are suitable over a large range of velocities. However, none of the 
prescriptions is universal. Depending on the nature of flow and geometry etc., one 
can always go for the best suited discretization procedure (see "Vanka, 1987, Fletcher, 
1988). It is also very difficult to be conclusive about any algorithm with its universal 
sense of applicability. In a recent investigation, Kim and Benson (1992) have compared 
three well known algorithms for modeling unsteady flows. They have concluded that 
a modified version of MAC algorithm (Harlow and Amsden, 1970; Braza et al. , 1986) 
is highly accurate and has some specific advantages on other schemes. However, the 
variants of MAC algorithm are not without problems of their own. Although, these 
variants have all usual disadvantages of the explicit methods, we prefer to use the MAC 
method because of the following advantages : it does not need boundary conditions for 
pressure, it can describe unsteady flows well, it is robust and uses primitive variable 
which are preferable in 3-D flow calculations. 


2.4 A Review on Turbulence Modelling 

i" i - ' 

2.4.1 General ’ ' 

Any flow whether laminar or turbulent, is fully represented by the Navier Stokes 
equations. The Navier Stokes equations can be solved on a fine enough grid with 
an exceptionally accurate discretization method so that both the fine scale and large 
scale aspects of turbulence can be calculated. This is termed as the Direct Numeri- 
cal Simulation (DNS) of turbulence (Rai and Moin, 1991). The modeling effort and 
simplifications that are employed in the study of turbulence are the consequences of 
the difficulty encountered in solving the full Navier Stokes equations on a fine grid. In 
numerical solutions, where DNS of turbulence is performed, the mesh spacings and the 
time steps need to be significantly less than those over which appreciable variation of 
velocity occurs, otherwise details of the evolution will not be correctly reproduced by 
the numerical prediction. The length-scale-range of the eddies of varying sizes (smallest 
eddies of the order of mm’s in the domain where the mean velocity is not much grt^ater 
than 100 m/s) and the time-scale-range of the velocity fluctuations due to the eddy- 
ing motion cannot be economically resolved by ordinary discretization methods. Tliis 
difficulty is circumvented through turbulence modeling. 

2.4.2 Various Approaches for ModeUing Turbulent Flows 

Launder and Spalding (1972) describe the modeling approaclu^s for turbuhuit flows 
comprehensively. The basis of turbulence modeling is to decompose the high fre(|uency 
fluctuating instantcineoiLS velocity of the fluid into a mean and a fluctuating ('ornpoiKuit, 
and solving for the mean velocities, while the effect of the fluctuating coinponcmts on 
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the mean motion is modeled after emidrical relations obtained for specific cases from 
experiments. The Navier Stok(5S equations together with the equation of mass conserva- 
tion form a closed system of ('(piations. Wiien the velocity com])onents are d('comi)osc'd 
we have more unknowns than number of available equations. The modified system of 
equations cannot be closed within itself unless empirical relations are supplied from 
experiments to correlate the fluctuating components with the mean motion. This is 
termed as the closure problem. 

The averaging of the equations to get tlie mean components can be done in various 
ways (see Hinze, 1987). The time averaging is suitable for the flows in which mean flow 
would be stationary or slowly varying in time. The ensemble averaging is basically a 
statistical concept wliich signifies tlie same meaning as time averaging for stationary 
flows, but is somewhat more meaningful for transient mean flows. A third alternative is 
the space averaging. The first two methods have been well tested and established while 
the tliird is used in a technique known as Large Eddy Simulation. The Large Eddy 
Simulation (LES) is another emerging area of tur bule nt research (see Schumann, 1975; 
Schmidt and Schumann, 1989; Robichaux, Tafti and Vanka, 1992 ). In LES, large scale . 
motions are calculated from DNS and subgrid scales are modeled. The most important 
aspect of this method is spatial filtering of the flow quantities i.e., decomposing the 
variables into resolvable and subgrid components before analysis. The interpretation of 
the results in physical terms is also quite tricky and complicated. 

The time averaged equations of motion completely resemble the original transport 
equations, except that they have a few additional terms containing the contribution of 
the fluctuating components of velocities to the transport mechanism. These terms are 
same for time average and ensemble average for stationary and homogeneous turbulence 
and the assmnption is known as ergodic hypothesis. Since turbulence is considered as 
eddying motion and the aforesaid additional terms are added to the viscous stresses 
due to mean motion in order to explain the complete stress field, these terms are called 
turbulence stresses or Reynolds stresses named after Osborne Reynolds. Space averaging 
culminates into a few more additional terms besides the Reynolds stress terms. 

Ferziger (1987) discusses in detail various alternative api:)roaches to turbulence mod- 
eling and explains the relative merits and demerits of these approaches. It is observed 
that many of the widely used models are based upon the eddy viscosity concept of 
Boussinesq. The most important effect of turbulence on mean flow is that it absorbs 
kinetic energy from the mean flow and increases the rate of transport of mass, momen- 
tum and energy in the direction normal to the streamlines of the flow. As both of these 
eflPects are mediated by the viscosity in laminar flows, it is natural to assume that the 
effects of tmbulence on the mean flow can be represented as an increased viscosity. This 
helps in expressing the turbulent stress as 
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where fit is called eddy viscosity and is attributed to random fluctuations of turbulent 
motion, k is the turbulent kinetic energy of turbulence and Sij is the Kronnecker delta. 

The eddy viscosity models are indeed popular, although they are criticized for assum- 
ing the turbulent eddy action to be isotropic by associating it to a viscosity. However, in 
order to determine the eddy viscosity, Prandtl proposed a mixing length hypothesis by 
analogy of eddy activity to molecular interaction (Schlichting, 1987). Thus, a need was 
felt for determination of a characteristic length and a velocity scale to define the eddy 
viscosity completely. The various models described in turbulent flow literature based 
on the eddy viscosity concept, necessarily attempt at determining a suitable velocity 
scale and a length scale to calculate 

The kinetic energy of fluctuating velocity components which defines the amomit of 
energy available for eddy action is proposed to be a good representative of turbulent 
velocity scale. The first step beyond the mixing length is such a model in which the 
transport equation for the above mentioned kinetic energy is derived. The length scale 
of eddies is continued to be derived algebraically. Such models are called as one equation 
models. 

The rate of dissipation of the turbulent kinetic energy (e) is subsequently proposed to 
represent the length scale of turbulence and the well known Prandtl-Kolmogorov relation 
is obtained as Ht = pC^k- fe. This gave rise to a family of two equation of models (Rodi 
and Spalding, 1970; Jones and Launder, 1972; Launder and Spalding, 1974 ) where 
the equation for turbulent kinetic energy determines the velocity scale. This equation 
can be derived from Navier Stokes equations by subtracting the mean equation from 
the unaveraged equations to obtain an equation for the fluctuating velocity. Taking 
the scaler product of this equation with the fluctuating velocity and averaging yields 
the equation for the turbulent kinetic energy. An equation for the dissipation can be 
derived from a Navier Stokes equation by an extension of the method used to derive the 
energy equation. The two equation models are quite successful and have become very 
popular. With only little modification, they are able to simulate a large variety of flows 
with reasonably good degree of accuracy. A comprehensive review of different kind of 
two equation models is available in open literature (Nallasamy, 1987). The equations 
derived in the above mentioned way apply to fully developed turbulence. 

The laminar region near a wall is handled by wall function treatment. In standard 
wall function approach (Launder and Spalding, 1974), the wall shear stress is related to 
the values of average velocity, kinetic energy and the distance outside the viscous layer 
using the generalized logarithmic law. Tliis wall function approach, by not requiring the 
integration up to the wall, luis a very beneficial effect on the economy of a calculation 
procedure. The viscous sublayer need not be integrated and therefore, a coarse grid can 
be used for calculation purpose, making it computationally economical. 

Although the thickness of the viscous zone is usually two or more orders of magnitude 
smaller than the overall width of the flow, its effects extend over tlie whok; How fi(!ld 
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since nearly 50 ijercent of the velocity change from the wall to the fix'estream occm’s in 
this region. Amano (1984) suggests a model where the near-wall region is divided into 
two distinct parts, a viscous sublayer region (0 < Y'* < 11) and an overlap layer (11 < 
< 400). In th(' viscous sublayer region, kinetic energy shows parabolic variation witli 
and in the overlap region it increases linearly. Amano (1984) proposes a tliree-layer 
model for the near-wall region that consists of a viscous sublayer (0 < Y''' < 5), a 
buffer layer (5 < Y"^ < 30) and an overlap layer (30 < Y'^ < 400). The behaviour of 
the turbulent kinetic energy, k and turbulent shear stress are proposed in the viscous 
sublayer and the buffer layer. Bakewell and Lumley (1967) suggest that within the 
viscous sublayer the streainwise fluctuating velocity increases linearly with the distance 
from the wall, y and increases approximately as in the buffer layer. The suggestions 
of Bakewell and Lumley are implemented in the near-wall three-layer model of Amano 
(1984). ■ 

Though the wall function approach is an economical method from the point of view 
of computer resources, the areas of steep gradients in the viscous sublayer are not 
resolved properly if the first point on the grid is placed away from the wall outside the 
viscous sublayer (Y^ > 11.63). The wall functions relating to the velocity and turbulent 
quantities at the first grid point depend heavily on the assumption of a logaritlimic 
velocity distribution and of local equilibrium of turbulence. These assumptions are not 
valid for separated flows and if there is a secondary flow extending into the sublayer. 
Therefore, considerable effort has gone into developing models that simulate the effect 
of viscosity and other turbulent quantities near the wall. Jones and Launder (1973) have 
extended the k — e model to the low Reynolds number form which allows a calculation 
right to the wall. Since viscous diffusion of k. and e are of the same order of magnitude as 
the tiubulence diffusion in the near wall region, the molecular viscosity is also included 
in the respective diffusion terms in different directions. A fmiction is introduced in 
the eddy viscosity equation to mimic the direct effect of molecular viscosity on the shear 
stress and the eddy viscosity is given by, fx-t = ffiC^pkr/e. The function ffj, is so chosen 
that it damps down from unity in the fully turbulent layer to zero at the wall. In all 
the models, is assumed to be a van Driest type of experimental damping function of 
local Reynolds number of turbulence {Rt = Ar/e Ut) and a similar independent variable 
{Rk = k^^~y/ut). Another fmiction /i is introduced in the production term of the k — e 
equations in order to increase the magnitude of e near the wall to account for the 
additional dissipation. In addition to these, a function /2 is used in the destruction 
term of e equation to indicate the low Reynolds number effects on the decay of isotropic 
turbulence. All these concepts are simplified and implemented in a somewhat more 
convenient way by Lam and Bremhorst (1981). From the evaluation of DNS data 
for channel flow a new model has been proposed by Rodi and Mansour (1993) which 
simulates the net effect of production terms and consequential increase in destruction 
term of the k — e model near the wall more accurately. 

An isotropic eddy viscosity model gives equal turbulent intensities in all directions. 
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The drawback of assuming isotropy is overcome by using Reynolds stre ss model. The 
Reynolds stresses are evaluated from exact transport equation for u'u'- which may be 
obtained by taking a velocity- weighted moment of the Navier Stokes equation and av- 
eraging. This will result in six partial differential equations. Many of the terms of this 
equation are not directly determinable and must be modeled. It is beyond the scope of 
the present study to give a comprehensive review on methods for closing the additional 
equations. But the contributions from Launder et al. (1975), Gibson and Launder 
(1978), Hogg and Leschziner (1989a, 1989b) indicate the successful implementation of 
this concept. However, it seenas that the Reynolds stress model is not without problems 
of its own. On the deficit side, Reynolds stress model continuous to be a computationally 
intensive option even on extremely powerful computers. 

An Algebraic Stress Model is proposed by Rodi (1976, 1984) which predicts anisotropy 
without requiring the solution of differential equations for these stress components. This 
approach begins with the exact transport equations for Reynolds stresses. Subsequently, 
it assumes isotropic dissipation, a modeled pressure strain term, and that transport of 
nX is proportional to the transport of k. Using these assumptions, the partial differen- 
tia/ equations for the Reynolds stresses are reduced to a system of algebraic equations 
at each point in the flow. Equations for the turbulent kinetic energy and dissipation 
rate are written. So, finally six algebraic equations are solved together with the k e 
equations. This model has considerable appeal and will keep on improving. It may be 
mentioned that the time scale for the Algebraic Stress Model is the same as other k-c 
models (k/e). 

The commonly used k — € models of turbulence are thought to be incapable of accu- 
rately predicting turbulent flows where the normal Reynolds stresses play an important 
role (mean velocity h as a gradient). Speziale (1987) has developed non-linear k-e model 
by adding a non-linear term to the linear Boussinesq approximation for the stress used 
in the standard k - e model. However, his non-linear model shows improvement both 
in the prediction of turbulent stresses and in the mean flow field. Speziale has also 
developed non-linear model equations to accoimt for the anisotropic behaviour of the 
' turbulent stresses without disturbing the tensorially invariant eddy viscosity structure of 
the two equation model. The model adds second order derivatives to the Boussinesq ap- 
‘ proximation of the turbulent stresses in the momentum equations. Dutta and Acharya 
(1993) extend Speziale’s non-linear model to heat transfer studies. The modified k — e 
model (Leschziner and Rodi, 1981) is another significant improvement on the standard 
k — e model by incorporating the effects of streamwise curvature. The improvement of 
the velocity related terms is expected to improve heat transfer. Therefore, it can be 
argued that the curvature effect should also be incorporated in the expression for the 
turbulent Prandtl number which in turn may bring about a change in the rate of the 
heat transfer. Dutta and Acharya show that the heat transfer prediction of the non- 
linear k — € model is found to be similar to that of the standard k — e model. However, 
heat transfer predictions of Dutta and Acharya with the modified k — ^ model agree 
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excrllcntly with oxiiorinKnital data of Aung and Goldstein (1972). 

As can be secui from the cited literature on turbulence modeling, no generalized 
turbulence model has yet been devised which can be apjdied through a modest com- 
putation effort. For turbulence models, comparison with experimental data is the only 
way to assess the accuracy and applicability of a model for a specific physical problem. 
Recently, Zhu, Fiebig and Mitra (1993) have computed a tliree dimensional turbulent 
flow with longitudinal vortices embedded in the boundary layer on a cliannel wall. They 
use standard k — e model for their simulation. Although the behavioiu oi k — e model 
is not known for such highly vortical flows, the comparison between the measurements 
of Pauley and Eaton (1988 a,b) and the computed results of Zhu, Fiebig and Mitra 
shows that the interaction of the longitudinal vortices with the boundary layer witliin 
a turbulent channel flow is captured quite correctly by the numerical simulation. 

It is envisaged that the A- — e model with some modifications in the wall function 
treatment will be able to predict our flow field which is also dictated by the longitudinal 
vortex induced swirling motion. For the validation of the model we also 2 >ropose to 
compare our results with the experimental data of Pauley and Eaton (1988 a,b). 


2.5 Remarks 

The present state of art of enhanced surfaces with respect to the augmentation of heat 
transfer has been outlined in the section 2.2 of this survey. The prediction of the details 
of the flow field over such enhanced surfaces and to determine the influence of different 
geometrical parameters on the flow and temperature fields by numerical solution of 
governing conservation equations are the essential goals of this study. Such prediction 
methods have been under development over more than two decades now. In section 2.3 
of the present chapter we have discussed the salient features of predictive procedures. 
The section 2.4 gives an account of various modeling techniques for resolving the issues 
related to turbulent flows. Apart from the references cited in the earlier sections, there 
is a host of other articles which will be discussed in the subsequent chapters together 
with the analyses presented therein. 



Chapter 3 

Mathematical Formulation 

3.1 Introduction 


Arising out of the practical and academic importance of the flow conflguration and 
thermal conditions as described in the previous chapters, a need is felt to compute three- 
dimensional laminar and turbulent flows in a horizontal channel with built-in wing-type 
obstacles. For detailed investigations, only one delta wing or a winglet-pair inside a 
horizontal channel has been considered (Fig. 3.1). The numerical simulation enables us 
to obtain a quantitative information on flow structure and heat transfer in the entire 
domain, with their dependence on the various governing parameters. 


3.2 Statement of the Problem 


As it has been stated above, computation is performed in a channel which is formed by 
two neighboring flns (Fig. 3.1). An obstacle in the form of a delta wing or winglet-pair 
(of zero thickness) is placed inside the channel. The base of the wing is flxed on the 
bottom wall and the apex faces the incoming flow stream at an angle of attack. In the 
case of winglet-pair, one side of each is fixed on the bottom wall and the trailing edge 
of each is free. Since symmetry prevails in the central vertical plane of the channel, the 
flow field in only half of the channel is computed. In the present study, the numerical 
simulations deal with two different regimes of flow, viz.^ laminar and turbulent. In the 
subsequent sections, the governing conservation equations for laminar and tvubulent 
flows and the appropriate boundary conditions will be discussed. 
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Type of obstacle ( protrusion ) 



Delta wing Winglet pair 


Figure 3.1; Flow model for computation 
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3.2.1 Governing Equations for Laminar Flow 

Tho dimensionless equations for continuity, momentmn and energy for laminar 
incompressible flow may be expressed in their conservativ(" forms as following ; 
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In the above equations, velocities have been nondimensionalized with the average incom- 
ing velocity Uav at the channel inlet, all lengths with the channel height H, the pressure 
with pU^.^, and the nondimensioiial temperature is defined as d = {T — roo)/(T^, — Too). 
The Reynolds number and Prandtl number are defined by Re and Pr respectively. 

3.2.2 Boundary Conditions for Laminar Flow 

The boundary conditions of interest in this investigation are (refer to Fig. 3.1): 


• Top and bottom plates: 

U = V = W = 0; 9 

• Side wall {z = B/2) and midplane ( z = 0 ): 


01-7 I 


1 
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9 = 0 

(3.8) 


w = 

At the channel inlet: 


• At the exit, a smooth transition through the outflow boundary is ensured by 
setting: 

dx^ ~ ax^ ax'- 'ax' ^ ’ 

No-slip boundary conditions for the velocities on the obstacles are used. The details of 
kinematic boundary condition on the surface of the vortex-generators will be discussed 
in a subsequent chapter. The temperature of the obstacle is considered constant and 
equal to T^^,. 
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3.2.3 Governing Equations for Turbulent Flow 

The time averaged Navier Stokes equations for turbulent incompressibl flow where 
the Reynolds stresses are expressed via the eddy-viscosity concept have been used as 
governing equations. These equations are written in a Cartesian tensor form as 
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where Uj and 9 are nondimensional time-mean velocity components and temperature, 
Ut^n is the nondimensional kinematic viscosity, is the nondimensional eddy diffusivity 
and kn is the nondimensional kinetic energy of turbulence. The subscripts i and j can 
take the values 1, 2 or 3 in the three coordinate directions Xi, X 2 and X 3 respectively. 

The ^" 1 , X 2 and Xz are equivalent to X, Y and Z in Cartesian coordinates. The 

Reynolds number Re is defined on the basis of average axial velocity at the inlet and 
the laminar viscosity as Re = UoHjv. The Prandtl number of the fluid is denoted by 
Pr. 

The turbulent kinematic viscosity Ut^n is given by 

= Cf,Rekl/en (3.13) 

In the above equation C^ is a constant (equal to 0.09) and is the dissipation rate of 
turbulent kinetic energy (normalized with respect to U^/H). The nondimensional eddy 
diffusivity at^,^ is given by 

= Cf^RePrkfJ {at e„) (3.14) 

where at is the turbulent Prandtl number. Generally, (j(=0.9. The modeled transport 
equations for k^ and are 
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Here, is the generation of turbulent kinetic energy (nondirmuLsioual) given by 

(1 + M/^^. dUXdUi 
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(3.15) 
(3. IS) 

(3.17) 


and ak,a^, Ci,C 2 are empirical constants : a,. = 1.0, < 7 , = 1.3, C’l = 1 . 44,63 = 1-92. 
These are based on wide range of experimental data. The argument for clu^osing these 
values are given in Launder and Spalding (1972). 
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Figure 3.2; Cross-section of the computational domain 

3.2.4 Boundary Conditions for Turbulent Flow 

The inlet conditions for velocity and temperature can be specified using the profiles 
we want to use. The turbulent kinetic energy kn and its dissipation rate are calculated 
from the value of turbulent intensity specified at the inlet. 

The inlet boundary conditions can be specified as following: 

U{Y) = {Ur,nHEY^))/X 
V =0 
W =0 
kn = 1.5 /“ 

^n{Y) = {kTC]J^)lxY for y < (A/x) 

= )/A for y < (A/x) 

where, Ur,n is nondimensional friction velocity, Y"^ is given by yur/i'^ I is turbulent 
intensity, x — 0.42 which is known as von Karman constant, A is a constant prescribing 
ramp distribution of miying length in boundary layers and equal to 0.09 and E=9.0. 

The boundary conditions for the outlet and the confining walls (both no-slip and 
free-slip) as shown in Fig. 3.2 are treated as follows: 
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• outlet: 


||:=0; f = {Uj,W,9,K,e,,) 
aX 

(3.19) 

• free-slip (symmetric) wall (plane I, III): 


W = 0, 11 = 0; / = (F, 7,0, 

(3.20) 

• free-slip wall (plane II): 


7 = 0, ^=0-, f = {U,W,9,K,en) 
oY 

(3.21) 

• no-slip wall (plane IV): 


U = V =W = 0,9 = 1 

(3.22) 


The wall functions due to Launder and Spalding (1974) and Amano (1984) are used 
to mimic the near wall region. 
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Instead of using equation (3.16) near tlie wall, the e„ at point p is computed from 
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In order to calculate the wall functions for temperature, the heat flux at the wall can 
be expressed in the following way; 


• For y/ <11.63 


• and for > 11.63 
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where, is the nondimensional heat flux at the wall given by 

Qw 
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It may be mentioned that Pfn in equation (3.31) is the Pee function which may 
be written as (see Dutta and Acharya, 1993) 
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3.3 Grid System Used 

Computational domain is divided into a set of rectangular cells (Fig. 3.3a) and a 
staggered grid arrangement is used such that the velocity components are defined at 
the center of the cell faces to which they , are normal (Fig. 3.3b). The pressure and 
temperatirre are defined at the center of the cell. In such an arrangement, pressure 
difference between two adjacent cells is the driving force for the velocity component 
located between the interface of these cells. The pressure field will accept a reasonable 
pres sure d istribution for a correct velocity field. Another important advantage of such 
a grid system is that transport rates across the faces of the control volumes can be 
computed without interpolation of velocity components. 
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Figure 3.3: (a) Grid spacing in the computational domain and the location of the 
Wing-type vortex grid showing the discretized variables (b) Three-dimensional staggered 
grid showing the locations .of the locations of the discretized variabk's 
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Figure 3.4: Boundary conditions and fictitious boundary cells 


3.4 Numerical Boundary Conditions 

So far, we have discussed boundary conditions for laminar and turbulent flows in 
symbolic form. In this section, we will briefly discuss the implementation technique of 
the boimdary conditions in the solution code. The boundary conditions are imposed 
by setting appropriate velocities in the fictitious cells surrounding the physical domain 
(Fig. 3.4). 


3.4.1 Boundary Conditions for Confining Walls 

Since, the equations are elliptic in space, we need boundary conditions on all confining 
walls, - even at the outlet. Consider, for example, the bottom boundary of the computing 
(physical) mesh. If this boundary is to be a rigid no-slip wall, the normal velocity on 
the wall must be zero and the tangential velocity components should also be zero. Here, 
we are considering a stationary wall. With reference to Fig. 3.4, we have 

= 0 

Wi^l,k = -Wi,2.fc 


for i 
and k 


2 to ire 
2 to kre 


(3.34) 


If the right side wall is a free-slip (vanishing shear) boundary, the normal velocity 
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must be zero and the tangential velocities should have no normal gradient. 


= 0 

UiJ,l = Ui,j,2 
= Vij,2 

The front plane is to be provided with inflow boundary condition. Any desired 
functional relationship may be recommended. Generally, normal velocity components 
are set to zero and a uniform or parabolic axial velocity may be deployed. Hence, with 
reference to Fig. 3.4, w« 

VFi,,-,, 

where, is the horizontal midplane. 

Continuative or outflow boundaries always pose a problem for low-speed calcula- 
tions, because whatever prescription is chosen it can affect the eiitire flow field upstream. 
What is needed is a prescription that permits fluid to flow out of the inesii with a mini- 
mum of upstream influence. In MAC family of algorithms, the second derivatives of the 
dependent variables in flow direction are set to zero in order to ensure smooth transition 
through the outflow boundary. These are implemented in the following way: 



= — Ui^ij'k 1 

for 

j 

= 2 to jre 



= 21/,;,. - 

> and 

k 

= 2 to kre 

(3.37) 


= 2W,j,k - Wi.,j,k J 


i 

= iim — 2 





3 can write 


-V2,> 

-W2J,k 

1.0 

1.5 [l - 


for 

and 


= 2 
= 2 


to jre 
to kre 


for i 
and j 


= 2 to ire 
= 2 to jre 


(3.35) 


3.4.2 Boundary Conditions for the Obstacle 

A special effort is required to satisfy the kinematic boundary conditions on the 
vortex generator. From Fig. 3.5, the angle of attack of the wing-type vortex generator 
can be expressed as /3 = tan \8Y/6X). Figure 3.5 further reveals that the U and V 
components of velocity fall directly on the surface of the obstacle in the region of the 
wing and are therefore set equal to zero. Implementation of no-slip condition for U 
component velocity at the side edge of the obstacle needs some manii)ulati<,)ns. From 
Fig. 3.6, it is found that Uij^k+\i 2 = 0. In other words, we can write 

(for i = ia to ib, j = 2 to ja and k = 2 to kc) 

It may be mentioned that the difference of count between kc and 2 is equal to 
that between ib and ia. The W component of velocity also needs manii)ulati()n for its 
zero value on the surface of the obstacle. From Fig. 3.7 


(3.38) 
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Wi.i.k 


Figure 3.5: Relative location of velocity components and wing on X-Y plane at Z—0 




Wij,k 

= -Wr-lJ,k 

‘■K 

f— 1 

1 

= -Wij,k 

Wi-,j,k 

1 

II 


(for i= ia to ib, j = 2 to ja and k = 2 to 


> 


kc) J 


(3.39) 


The temperature boundary conditions on the wing-type obstacle are diagram- 
matically shown in Fig. 3.8. The obstacle boundary conditions for temperature are : 


di,j,k = l,j,k 

^ij,k = 20 ^, — dij-l,k 

^i-l,j+l,k ~ 29^^ lj,fc 

~ 29t^ 9i^2,i+\,k 6tC. 

(for i= ia to ib, j = 2 to ja and k = 2 to kc) 


(3.40) 


In practice, the vortex generators can easily be manufactured by punching or 
embossing the wall. In some computational results, the effects due to the hole beneath 
the vortex generator have been taken into account. It is considered that the heat 
exchanger cores have a number of plates and on each plate the vortex generators are 
punched out in such a way that the holes are perfectly aligned. In our computational 
domain, a small part of top and bottom boundaries must be set to identify the punched 
holes. For this purpose spacewise-periodic boundary conditions on the location of the 
punched holes have to be applied in the y - direction (Fig. 3.9). Here, the period length 
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Y-Z plane 




Figure 3.6: 


Side edge of the wing-tyjie vortex geuerutor on Y-Z plane 
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j = jim 
j = jre 

j = ja 


j = 2 

j=l 


Dn of the spaiiwise velocity components and the wing-type 
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Figure 3.9: Periodic boundary conditions for pressure and velocities on the top and the 
bottom wall 


is channel height H. However, the boundary conditions for the fictitious cells on the 
bottom are : 


and on the top 




Vi,u = 

0-5(Vijre^k + 

= 


— 

Pi,jre,k 

(for all 

i,k in the punched holes) 


U'iyjim^k 

= 

Vijira,k 

= ViOk 

V • L 
^ re, k 

= 0.5{Vijre,k + 1^,1, fc) ^ 



P i,jim,k 

= Pi,2,k 

(for all 

i,k in the punched holes) ^ 


(3,41) 


(3.42) 


Such boundary conditions have also been discussed by Hirt et al. (1975) and 
Biswas and Chattopadhyay (1992). However, for plate-fin heat exchangers, a dif- 
ferent manufacturing technique for mounting the vortex generators is adapted. 
For this case, no hole exists underneath the vortex- generators in order to avoid 
mixing of the hot and cold streams. 
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3.5 Method of Solution 

A modified version of Marker and Cell (MAC) method due to Harlow and Welch (1965) 
and Hirt and Cook (1972) is used to obtain the numerical solution of the governing 
equations. 

3.5.1 MAC Algorithm 

Because of staggered grid arrangement (Fig 3.3(b)), the velocities are not defined 
at the nodal points. Whenever required, the velocities at the nodal points are to be 
found by interpolation. For example, with uniform grids, we can write Ui\i/ 2 ,j,k = 
0.5(17i.i ij,jt + Uij^k)- Where a product or square of such quantity appears, it is to be 
averaged first and then the product is to be formed. 

Convective terms of the momentum equations are discretized using a weighted 
average of second ‘upwind and space centered scheme (Hirt et al, 1975). Diffusive terms 
are discretized by a central difference scheme. Let us consider the discretized terms of 
aj-momentum equation (refer to Figure 3.3(b)). 


(3.43) 


(3.44) 


(3.45) 


dU- 


dX 




4SX 


[{Ui,j,k + LVij./t) (Uij,k + Ui+ij^k) 


Oip I (^LiJ^k ” 1 ” L i-i-lj^k) i {Lij,k Ui+lJ^k'} 
— {Ui-ij^k + Ui.j,k){Ui^ij^k + Uij^k) 

Up I (Ui-ij^k 4“ Cij^k) I (Ui^ij^k Hij^k)] 

= DUUDX 


dUV 


dY 


■J 


A6Y 


[(Lij,fc + Vi+l,j,k) {Uij^k + Uij+l^k) 


+ Up I I {U,,k - ^ 7 .;>U-) 

~ l.k){Uij-l^k p 

- Up I I {UiJ.Uf: - U,J,k)] 

= DUVDY 


dUW 

dz 


J ij,k 


46Z 




+ Up I (Wij.fc + W^Hju) I {Uij,k - 

- {^i,i,k-{ + ^ik\,j.k-\)[Uij^k-\ 4- L^,,;,fc) 

“ Up I (Wij't.i -h Wnij^k-i) I {Ui,j,k-\ - ih.j,k)] 

= DUWDZ 
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dP 

dX 

d'Tr 

dX^ 

d'^U 

dKr 

dz^~ 


= DPDX 


^ ’ i 4 3 ,j,A: 2 Ui - ] J^k 

(6Xj^ 

^^2.j4-l,fc ^^h,j,k UjJ -l,k 

(6y)2 

Z , t = D2UDZ2 


= D2UDX2 


= D2UDY2 


(3.4(5) 

(3.47) 

(3.48) 

(3.49) 


with, Op — > 1, the scheme tends to Second upwind and Op — ^ 0, the scheme tends to 
Space centered. 

The factor Qp is ciioscn in such a way that the differencing scheme retains ’something’ 
of second order accuracy and the required upwinding is done for the sake of stability. 
A typical value of ap is between 0.2 and 0.3. 

Solutions for the velocities are obtained in two folds. First the velocity components 
are advanced explicitly using the previous state of the flow, having calculated acceler- 
ations caused by convection, diffusion and pressure gradients through a time step 6r. 
The explicit time increment may not necessarily lead to a velocity field with zero mass 
divergence in each cell. In the subsequent fold, adjustment of pressm-e and velocity is 
done by an iterative process in order to ensure mass conservation in each cell. If the 
explicitly advanced provisional velocity is termed as Ul^l, then from equation (3.2) we 
can write 

C’A; = urj.t + StIRESIDU]1,_i, (3,50) 


where 


[RESIDUyL^ = [{-DUUDX - DUVDY - DUWDZ) - DPDX 

+ {l/Re){D2UDX2 + D2UDY2 + D2UDZ2)] (3.51) 


In a similar manner we evaluate 




Km + St[residv] 


n 


(3.52) 


= IKm + ^r[RESIDWyL, (3.53) 

As discussed earlier, the exphcitly advanced tilde velocities may not necessarily 
lead to a flow field with zero mass divergence in each cell. This implies that at this stage 
the pressure distribution is not correct. Pressure in each cell will be corrected in such a 
way that there is no net mass flow in or out of the cell. In the original MAC method, the 
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corrected pressures are obtained from the solution of a Poisson (Kpiafion for pressure. 
A related technique developed by Chorin(19()7) involved a siniultau<>ous iteration on 
pressures and velocity components. Vicelli(1971) has shown that tlu'. two methods 
as applied to the MAC method are equivalent. We shall make use of the iterative 
correction procedure of Chorin (1967) in order to obtain a divergence-free velocity field. 
The mathematical methodology of this iterative pressure-velocity correction procedure 
will be discussed herein. 

The relationship between the explicitly advanced velocity component and the veloc- 
ity at the previous time step may be written as 






+ 


- p’\. 


iM 


6X 


+ St[CONDIFI%^,, 


(3.54) 


[CONDIFUrij^, = [RESIDUAL, -f- [DPDXyL, (3.55) 

where, on the other hand, the corrected velocity-component (unknown) will be related 
to the correct pressure (also unknown) in the following way: 




Vtj.,+ST 


[0714-1 _ pn4«l 
l'^ iyj,k ^ if Ij,A: 


6X 


+ 6T[C0NDIFl%_f, 


From equations (3.55) and (3.56) we can write 


(3.56) 


Uij,, - U,j, - ^ 


where the pressure correction may be defined as 


(3.57) 


P' _ p»ifi 


pn 


(3.58) 


Neither the pressure corrections nor are known explicitly at this stage so one 

can not be calculated with the help of the other. Calculations are done in an iterative 
cycle and we write 




^r[PU 




sx 


In a similar way, we can formulate the following arr 9 


■V: 


(3.59) 


Corrected 

Estimated Correction 



sx 

(3.60) 

rjn+l 


sx 

(3.61) 

- 

hj,k 


+ tBiAzIUiiA 

f)Y 

(3.(i2) 
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T/7.+ ] 




ij-hk 

Wl I MPUk - PjiMi] 

'^i.i.k err 


W’^D W“+^ _ ^i,j.k-\. 

” i,j,k-l cr^ 


(3.63) 

(3.64) 

(3.65) 


The correction is done through continuity equation. Plugging-in the above rela- 
tionship into the continmty equation (3.1) we shall obtain 

rc/-;{- erv,-., wrf_i-wsi.A 

6X 6Y 6Z 

6X 6Y ^ 6Z 

. \ Pl».u-‘^Plu + PUi.k . PlMi-iPU + PlMJ. ' 

I (ixy- (SYy 

P' —OP' 4 - P' 

+ (3.66) 


usi - u;^,.k - ws;., 

«x iy 6Z 

u^tk - , Ki* - 

6X 6Y ^ 6Z 


\ 26t{PI,,) 26t{PL,) 2St{P!,^,: 

^ 6X^ ^ <5y2 ^ 6Z^ 


(3.67) 

(3.68) 


In deriving the above expression, it is assumed that the pressure corrections in 
the neighboring cells are zero. Back to the calculations, we can write 


0 = (£))ij. t + P* 2STf— + — + — 


SZ^/J 

+ (5,69) 

+ 5^2 + iZ-VJ 

In order to accelerate the calculation, the pressure correction equation is modified 


_i_ + J_ 
5^2 ^ sz^ 


(3.69) 


D/ ^o{D)ij,k 

o c ( 1 I 1 I 1 

-51^ >5^^ 


(3.70) 
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where cJo is the overrelaxation factor. A value of a;„ = 1.7 is coiumouly used. After 
calculating the pressure in the cell (i,j,k) is adjustc'd as 


pn+l 




(3.71) 


Now the pressure and velocity components for each cell are corrected through an 
iterative procedure in such a way that for the final pressure field, the velocity divergence 
in each cell vanishes. The process is continued till divergence-free velocity field is reached 
with a prescribed upper bound; here a value of 0.0001 is recoiimiended. This solution 
scheme is continued until a steady fiow is obtained. 

Finally we discuss another important observation. If th<‘ vt'kxdty boundary con- 
ditions are correct and a divergence- free converged velocity fi<'ld has Ixm'u ol>tained, 
eventually correct pressure will be determined in all the cdls iuchxiiug tlu^ cells at the 
boundary. Thus, this method avoids the application of pressun* boundary conditions. 
This typical feature of modified MAC method has been discuss(Hl in inoixi d<?tails by 
Peyret and Taylor (1983). 


3.5.2 Discretization Technique for k and € Equations 


Convective terms of the k and e equations rn'e also discretized using a wc'ighted average 
of second upwind and space centered scheme (Hirt et «/., 1973). Diffusive tc'rms are 
discretized by a central difference scheme. 

Let us consider the discretized terms of the k equation (nder to Figure 3.3(b); here 
the location of P or 0 is similar to the location of k) 


dkU 

dX 




26X 


+ kn-ij^k) Uij^k 


+ ocp{ki,j,k - A’i+i.i.fc) I Ui,j,k I 

~ (^'t F kij^k)Ui~]j^k 

- rtp I I {ki-ij,k ~ ^■ij./t)] 

= DkUDX 


(3.72) 


■ d 


dX 

\ 'dXJ^ 


2(6X)- l.j.A:) (^'il \.j,k ■" ki^j^k) 

i.j,k + {i^t)ij.k){k,j,k - k , i,/,t)] 

= DIFFkX 


(3.73) 


The difference (luotients involved in the Generation term aio disciotized as following. 

9 IL — - Uj i.i.k nr.-m.- /-> 'tO 

OX 6X =^UDX (3.74) 
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dV 

dV 

dW 

W 

dU 

dY 


dV 

dX 


dW 

dX 


dU 

dZ 


dV 

dZ 


dW 

dY 


= ]^VDY 

= DWDZ 
bZ 


W^i,hk + Ui-lj^k + Uij+\^k + Li^ij+I^k) 


AbY 


{Ui,j,k + Ui--ij^k + Uij-]^k + 

DUDY 

+ Yi+\,j,k + 

+ Yi-l.j.k + 

DVDX 

7^ W^.j.k + + Wi+i,,-, + 

4oA 

+ W,j^k-1 + ■,_!)] 

DWDX 

[i^hj,k + + Uij^k+] + Ui^ij^k+l) 

{Ui,j,k + Lfi-lJ^k + Uij,k-1 + ^7i-lj,A:-l)] 

DUDZ 

[iYj,k + Vi,j-hk + ^ij,fc+l + Vij-hk+i) 
{Vij,k + ^i,j-l,k + ^i,j.k-l + K;-!,*-!)] 

DVDZ 

-^Y ^ij,k-l + Wij+I,k + Wij+i^k-l) 

{Wij^k + Wij^k-l + Wij-i^k + 

DWDY 


(3.75) 

(3.7G) 

(3.77) 

(3.78) 

(3.79) 

(3.80) 

(3.81) 

(3.82) 


In terms of the difference quotients, the e equation consists of the similar quotients 
as the k equation. 


3.5.3 Solution of the Energy Equation 

After evaluating the correct velocities, the energy equation is solved with an Successive 
Over-Relaxation technique to determine the temperature field. We shall discuss the 
solution procedure of the energy equation in this section. 

The steady state energy equation, neglecting the dissipation term, may be written 
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in the following conservative form as 

due dve owe _ v-e 

dY ^ dZ ~ Pe 


(3.83) 


Equation (3.83) may be rewritten as 

= Pe[CONVT]^jk (3.84) 

where, [CONVT]ij^k is the discretized convective terms on the left hand sid(' of ('(^nation 

(3.83) and m stands for iterative counter. To start with, we can as.sume any guess value 
of 9 throughout the flow field. U, V,W are known from the solution of luonumtum 
equation hence equation (3.83) is now a linear equation. Howevcu*, from tin; guess 
values of 9 and known correct values of U, V, and W the left hand side of equation 

(3.84) is evaluated. A weighted average scheme may be adapted for discretization of the 
convective terms. After discretizing and evaluating right hand side of equation (3.84) 
we obtain a Poisson equation for temperature with a source term on tlw^ riglit hand side. 
Now, we can follow a SOR technique for solving equation (3.84). Consider a discn^tized 
equation as 

9i+ij,k ^e^jk ~i~ 9{j^\^k ‘^eijk ~i~ e^ j^ik ~~ ^e^jk t e^ jk-i ntm 

{6xy- (dy )2 {6zy^ 


where = Pe[CONVT]”^jj^ 
or, 


or. 


2 2 2 

[(5X)2 ^ (6y)2 {6zy- 

= 


2 

J6xy + (dzy 


(3.85) 


where. 


ei+[,j,k + + 9ij . 1^^, Oij^ki-i + eij^k I 

{ 6 xr- (6y)2 + (^szy 


9ij,k in equation (3.85) may be assumed to be the most recent value and it may Ix' writt(ui 

bi order to accelerate the speed of computation we introduc(^ a ov<'rrelaxation 
factor w and 

= (3.86) 

where ^ is the previous value, is the most r(x:ent valium and B’"\} is th(‘ ealculat<‘d 
better guess. The procedure will continue till the required coimu'ge'ucT is achievcui. 
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For accuracy, the mesh size' must be cliosoii small enough to rt'solvc' the cxijected 
spatial variations in all dependent variables. 

Once a mesh has been chosen, the choice of the time increment is gov(>rned by two 
restrictions, namely, tlie Courant-Friedrichs-Lewy (CFL) condition and the restriction 
on the basis of grid-Fourier numbers. 

According to the CFL condition, material cannot move through more than one cell 
in one time step. Therefore, the time increment must satisfy the inequality 


6r < min 


f 6X 

6Y 

6Z 1 


V|’ 

|W J 


(3.87) 


where the minirnmn is witli respect to every cell in the mesh. Typically, 6t is chosen 
equal to one-fourth to one-third of the minimum cell transit time. When the viscous dif- 
fusion terms are more important, the condition necessary to ensure stability is dictated 
by the restriction on the grid Fourier numbers, which results in 


1 (sxr-jSYnszr 

2 {6X^- + + 8Z^) 


(3.88) 


The final 8t for each time increment is the minimum of the 6t’s obtained from equations 
(3.87) and (3.88). 


3.6 More on Discretization of Governing Equations 

We have already seen that the continuity equation is discretized by a simple first- 
order differencing in such a way that the mass conservation is maintained in each cell. 
The temporal terms and the terms corresponding to the pressure gradients are also 
discretized by first order differencing. The discretized viscous terms are second order 
accurate. All these terms have been presented in the previous section. The convective 
terms have been discretized by using a weighted average scheme wfiich combines upwind 
and central diflferencing to achieve the stability of the upwind method and better formal 
accuracy of the central differencing. We have already discussed that the factor, Op serves 
to balance the upwind contribution. If 0, the difference equations are centered in 
space. In the present investigation, the factor etp is restricted between 0.2 and 0.3 so that 
the formulation can retain "something" of the second order accuracy (Roache, 1972). 
Rimchal and Wolfstein (1969) and Timin and Esmail (1983) used this discretization 
scheme to compute driven cavity flows. Their results compare favorably with those of 
the second order schemes. 
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Vjjtl.k Vj+ij+i^k 
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Ui-I, 
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^i.j.l.k 

'ji.i.h 

Ui,j + I,k 
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^^l,j + l,k 
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«l>H,j,k 
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Figure 3.10; Dependent variables on a rectangular mesh for describing QUICK scheme 

3.6.1 Some Improvements on Discretization of Convective Terms 

Udimtidnluh ic? T' MAC method with a mul- 

alvoriZ for tht lm tliat their marching 

term d(U,l,)/aX may be represented’^ ‘hff"s<'d. Tlie convective 




(8.89) 
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where, the variables defined as 


^ j,A: 


2 l(e«, .,.1 + (['«.*, 


and 


(U^h 


■ i J'>c 


[{U(f>)i-ij,k - 2{U(j))ij^k + for Uij_k > 0 


(3.90) 


3 

1 
2 

“ - 2{U(l))i^ij^k + {U^)ij^k] for Uij^k>0 (3.91) 

The parameter ^ can be chosen to increase accuracy or to alter the diffusion-like char- 
acteristics. It may be pointed out that ^ = 3/8 corresponds to the QUICK scheme of 
Leonard (1979). However, the discretized second term of the .'c-momentum equation 
after the QUICK sclieme may be written as 

'd{uvy 


dY 


0.25 

86Y 


[3f/jj4.i^t(l/,7+l,t d" ^i,j,k d" U+lj'+l.r- d“ U+lj,fc) 


d- TUij^k{Yij,k d- Vi,j-l.k + "K+lJ./t d- U+l,;-i,i-) 

“ d- Vi,j-2,k d- d- K+l,j- 2 ,(t) 

d- Uij.2,kiYij-2,k d- Vi^j_2,k d- Vi+l,;-2,i d" U+lj-3,A:)] 

for Uij^k > 0 


(3.92) 


and 


dUV 


dY 


0.25 

86y 


— Uij+2,k{Vi,j+2,k + y,;+l,fc d- Ui+ij+ 2 ,fc d- 1 / 4.1 j+i,*) 


d- 7Uij+l^k{Yi,j+l,k + ^i,j,k d- Vj+lj’+l,* d- I/+l,_7,t) 

8Uij^ki.^i,j,k d“ d“ I/+l,j,t d“ I/+l,_7- l,t) 

— SUiJ-l^k{Yi,j-l,k + yi,j-2,k d- I/+l,;-l,fc d- Ui+lj_2,fc)] 

\ for Uij,A:<0 .r...,, ’ (3.93) 

We have combined the QUICK scheme of Leonard and a upwind scheme in order 
to improve the accuracy of discretized convective terms of the governing equations. 
The upwind discretization of the second term of the ^-momentum equation has been 
accomplished in the following way: 

'dUV 


dY 


up 


+ I {Vijj , + I (Ui^t - 

— (Kj-i.t + + Uij^) 


(3.94) 
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Finally, the discretization of the term d{uv)/dy for the .x-momeiitum equation has 
been accomplished in the following way: 


d{UV) 


dY 


= (1 - Oiy) 


d{UV) 


dY 


+ a. 


quick 


d{UV) 

dY 


(3.95) 


up 


where, ay is contribution of the upwind factor.The factor ay can be found in the followinff 
way: ' ° 


dVii 


Re% - 
Re’l'^ - 


(3.96) 


II 

fl — 1, 

1 KK 1 

1 W6Z r 


M (/ M 

r - 

1 // 1^ 

1 i> \ 


'l IfhX 1 

1 ^ 1 

1 WhZ |1 

J 1/ n 

1 1’ 

1 IJ 


where, 

(minimum cell Reynolds number) = min 

^gmax (maximum cell Reynolds niunber) = mfi:. 

Re\ = calculated cell Reynolds number of any cell. 

n ^ third-order upwind scheme is represented by blending the fourth- 

order central difference of the first derivative and the fourth-order dissipation term. As 
n examp e, he discretized second term of the .'c-momentum equation becoiru^s 
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3.6.2 Comparison of Results based on Model Problem 

Tlircc different discretization schemes wliicli iiavc bt'cn already discussed in the 
earlier sections are tested on a model problem of flow in a lid-driven square cavity. In 
this section, the results obtained by the above mentioned discretization schemes are 
reported and comparisons are made with the available results. 

Tlie computational domain consists of a two dimensional lid-driven square cavity 
with no-slip and impervious boundary conditions at the bottom and side walls except at 
the top, where u is a non-zero constant and equal to ULID. In the governing equations, 
the velocities have been nondimensionalized with respect to ULID. All lengths have 
been scaled with respect to cavity height H, and the Reynolds number is defined as 
Re = (ULID. iir)/i/. In order to display quantitative information on the velocity field, 
the velocity vectors for the above mentioned flow field for tliree diflferent discretization 
schemes have been shown in Figs. 3.11, 3.12 and 3.13 

In Fig. 3.14 u velocity components at the vertical mid plane of the square cavity for 
a Reynolds number of 400 are plotted. A (42 x 42) grid is used and the results due to 
three different discretization schemes are compared with the results obtained by Ghia 
et al. (1982). It may be mentioned that Gliia et al. (1982) have employed a 129 x 129 
grid via a multigrid technique. However, Fig. 3.14 shows reasonably good agreement 
between Ghia et aVs result and all the three cases of present computation. Even 
otherwise, the positions and the values of the minimum u velocity in position and value 
of all the computations agree with each other on an overall assessment. The v velocity 
components at the horizontal mid plane for the same Reynolds number (i?e=400) are 
shown in Fig. 3.15. It is evident that the results obtained due to all the three schemes 
demonstrate reasonably good agreement with the results due to Ghia et al. (1982). 

This has also been observed that the weighted average scheme takes considerably 
less CPU time to compute (discussion on computational time has been presented in a 
subsequent section). The computational accuracy of the weighted average scheme is 
also calibrated by comparing the peripheral average local Nusselt numbers predicted 
by our numerical scheme with those of available literature (Shah and London, 1978) 
for thermally developing and hydrodynamically developed flow through parallel plates. 
The results are found to be within about 3 percent for a Reynolds number of 500 and 
uniform wall temperature boundary condition. As such, in the present investigation, 
we have used weighted average schemes of discretization for computing laminar flows. 
The turbulent flows have been computed by making use of QUICK-upwinding scheme. 
For laminar flows, weighted average scheme has been found to be time optimized and 
reasonably accurate. 

We may mention that for laminar flows we have also compared the results due to 
our weighted average scheme of discretization with those due to QUICK discretization 
scheme for the problem of our investigation. The combined-spanwise-average Nusselt 
number (Nusa) distribution and the combined-spanwise-average skin friction coefficient 
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Figure 3.12: Velocity vectors in a lid-driven square cavity for i2e=400 (42 x 42 grid); 
discretization due to QUICK scheme 
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Figure 3.13: Velocity vectors in a lid-driveii simare cavity for, Re==m) (42 x 42 grid) 
discretization due to third-order npwiudiag scheme 
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Figure 3.14: Variation of U velocity along the vertical midplane for the lid-driven flow- 
in a square cavity 
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Figure 3.15: Variation of V velocity along the horizontal midplaue for the lid-driven 
flow ill a square cavity 
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Figure 3.16: Comparison of weighted average scheme with QUICK 

{Cf X Re) distribution in a rectangular channel with a built-in delta wing are the 
parameters which form basis of this comparison. The comparison is shown in Fig. 3.16. 
Maximum discrepancy in the prediction of Nuga is found to be less than 3 percent. 


3.7 Spatial Grid Independence 

An effort is undertaken to obtain grid independent results. For Re= 500 and Pr = 0.7, 
in a channel (a = 2) with a built-in delta wing (A = 1) at an angle of attack {j3) of 26'^, 
reasonable agreement between the results obtained for 15 x 26 and 20 x 30 cross-stream 
grids is found. In the foregoing specification of cross-stream grids, 15 and 20 refer to 
the number of gi'ids in the y direction. Similarly, 26 and 30 refer to the number of grids 
in the direction. It may be mentioned that for a J.5 x 26 cross-stream grids, 60 grid 
nodes are taken in x direction for a nondimensionaf length of 8.4. 
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Computations are also performed for a shorter channel of nondirrumsional length 
4.3 for three different grids. Other geometrical and flow paranuiters concerning these 
computations remain same as the earlier one. Computations for grids 30 x 1() x 13, 
30 X 15 X 26 and 30 x 30 x 52 have been carried out. The grids are specified in x, y and 
z directions respectively. However, the extrapolated grid-independent average Nusselt 
number differs from the corresponding Nusselt number obtained with the coarsest grid 
by less than 10 percent. A grid size of 30 x 15 x 26 produces average Nusselt number which 
differs from this extrapolated grid-independent average Nusselt number by less than 3 
percent. Hence most of the computations are done by using a 15 x 26 cn^ss-stream grids. 
In some cases, cross-stream grids finer than (15 x 26) are used for spc'.cial requirements, 
viz., accommodating some typical obstacle position and obtaining a clos('r point near the 
wall in the case of turbulent flow computations. Tlu; grids in tht^ str<'amwi.se directions 
are suitably chosen to cover the required length of the chaniufl with a reasonably good 
accuracy. 


3.8 Computational Time 

In the problem of a channel with a built-in delta wing, for a grid size of 60 x 15 x 26 
(in X, y and z directions respectively) and Reynolds number of 1000 (laminar flow) the 
CPU time with respect to steady solution is 157 min. 

For a Reynolds number of 15000 (turbulent flow), the corresponding CPU time with 
respect to steady solution is 620 min and the grid size is 30 x 15 x 27. 

The laminar flow computations are performed on a HP-9000/850 serif's computer 
whereas the computations for the turbulent flows are performed on a CONVEX C-220 

computer. The CONVEX C-220 consists of 2 processors with 128 MB RAM and 2 GB 
disk space. 



Chapter 4 


Evaluation of Performance 
Parameters for Laminar Flow 

4.1 Introduction 

Transport-enhancement in fin-tube and plate-fin heat exchangers is the essential goal 
of this investigation. Improvement in heat transfer is attributed to the decrease in 
thermal resistance owing to the spiraling motion of the flowing fluid. Thus the main 
quantities of interest are tlie enhancement of heat transfer and the corresponding flow 
losses. In this study a detailed numerical simulation has been done on an element of 
an heat exchanger which is representative of a part in the core region of the above 
mentioned heat exchangers. This chapter presents the results of the simulation for flow 
and heat transfer in a channel with both the vortex generators, viz., delta wing and 
delta winglet-pair. The first part of this chapter presents the performance of a delta 
wing as a vortex generator. In order to evaluate performance of the winglet-pair as the 
promoters of enhanced heat transfer, the latter part of this chapter presents simulation 
of flow and heat transfer in the configuration of a channel with built-in winglet-pair. 
The results presented in tliis chapter are in the laminar-flow regime. 


4.2 Performance of Delta Wing- Type Vortex 
Generators 

A number of computations have been performed at different angles of attack (/3) with 
a delta wing as a vortex generator. For slender delta wing in an infinite medium, leading 
edge vortices have a spiraling structure and they diverge slightly. Now, when we envision 
a wing moving in an infinite medium, the wake grows longer and basically becomes a 
swirling flow supported by trailing vortices. However the flow in a channel with an 
attached delta wing shows no trailing edge vortices. Figure 4.1 shows the generation of 
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vortices and their deformation as they move along the channel in ]>r('senc(^ of a delta 
wing. Deformation takes place due to the reduction in strcuigth of tlu^ vortices which is 
brought about by the viscous resistance of the wall. Tlie cross stream vcdocity vectors 
at different axial locations are shown in this figure. 

Figure 4.2 illustrates the static pressure distribution on the cross planes at the same 
axial locations for the same geometrical and flow parameters as those of Fig. 4.1. The 
nondimensional static pressure (p/pU;^) contours follow the same qualitative trends as 
those observed in experiments conducted by Hummel (1978). 

Development of streamwise vorticity, normalized with the mean axial velocity and 
the channel height, is presented in Fig. 4.3. Dotted lines indicate' ne'gative vorticity. 
Peak vorticity occurs in the vortex centers. However, the' pe'ak vorticity in th(^ vortex 
core decreases in the downstream direction. Reynolds numlxu- for tins case^ is 1000 and 
the angle of attack of the wing is 20'^'. 

The isolines for the axial velocity at two different axial locations bediind the wing are 
shown in Fig. 4.4. On the base plate behind the wing, the vortice's leuid to a thinning 
of velocity boundary layer in the downwash region at the; wing symme'try plane and 
thickening at the array symmetry plane. 

Figure 4.5 shows the isotherms over eight cross-planes locat('d at (ught diff(U'ent axial 
locations in a chamiel for Re =1500 and Pr = 0.7. It is sc'en that as the fluid stream 
moves from an axial distance X = 1.268 to X = 5.706, the vahu' of the isotherms in 
the core region increases and consequently the bulk ternpcu’aturc' k('('ps on rising. The 
relative location of the wing has been shown clearly in this figuri'. Due to the spiraling 
structure of the flow behind the wing, there is a mixing of cooh'r stri'am of the core with 

ttius, as can be observed^ the core t:(un|)(>ratiire rises in 

the downstream. 


However, the heat transfer and the skin friction coefficients, the major performance 
parameters are dependent on the type of the obstacle (wing or winglet pair), the 
eynolds number and the pometrical parameters of the obstacle. In the following 
sec ions we s all discuss the influence of these governing iiaranu'ters on the performance 

assumed as the working fluid; lu'iice the Praudtl 
number of this study is fixed at a constant value of 0.7. 

distiucti,,.! „( tl.,. Haa.fdr ,,«formau«^ the 

combined spanwise Nusselt number 


Nttsa — • 


2 Jo ^'iTwi{x,z)-Th{x))dz + 2jf\T,,,2(^^^^ - Th{x))(b 




position'^inthrdmnntl.'''^'^^ longitudinal location to depict heat transfbr at any axia 
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Figure 4.1: Cross-stream velocity vectors at different axial locations behind the wing 
showing generation and deformation of vortices in the channel 
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Figure 4.3: Contours of normalized streamwise vorticity at different cross-sections of 
the channel in presence of a delta wing 
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Figure 4.4: Isoliiies for axial velocity at two different (I’oss-stn'am plaiu^s (a) 0-58 
wing-chord downstrecun and (b) 1.74 wing-chord downstn'ain of a d<‘lta wing 
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Figure 4.5: Isotherms at different cross-planes in the channel in presence of built-in 
delta wing 






66 


CHAPTER 4. EVALUATION OF PERFORMANCE 


with higher strength which results in better heat transfcir. At a uoiuiiinensional axial 
location of 4 from the inlet, for an angle of attack of 25'' , observe an enhancement 
of 37 percent in combined spanwise average Nusselt number ovc'r tlu' case of channel 
without obstacle. Furthermore, at the same location, for /i=35, an improveirient of 8.87 
percent in Nusa is obtained over the case with l3 2o . 

In the above mentioned case, a fully developed velocity profile was considered at the 
inlet of the channel. For a developing flow, the effect of flow development in the channel 
and the effect of vortex generators on the flow field are superimposed in the downstream. 
In the earher case, in order to distinguish the effect of vortex generators (to eliminate the 
effect of flow development), a fully developed velocity profile at the inlet was deployed. 
However, here we shall include the effect of developing flow and the influence of punched 
holes beneath the vortex generators on the heat transfer performance. 

Figure 4.7 shows the comparison of cross-stnuun velocity v('ctors at diih'n'.nt axial 
locations from the inlet of the channel for the cases without and with stamping on 
the solid walls. Due to the stamping on the walls, a velocity field normal to the vortex 
motion is induced in the downstream direction. This induc(Hi downward normal velocity 
field reduces the strength of the vortex flow and as compHnnl with th<' case where there 
is no stamping, a decayed circulatory flow patt(U’n is obs(;rv('d at th<' sanu' axial location. 

Figure 4.8 shows the heat transfer performance in a channel for a simultaneously 
developing flow (Kakac, Shah and Aung, 1987) in presence of a d('lta-wing. Figure 4.8 
also compares the heat transfer on the channel walls for the with and without stamping 
cases. For the case of without stamping, in the region of the wing (from X — 2.92 to X = 
3.79), the combined spanwise average Nusselt number rist's to a lugh valium of 10 and then 
takes a plunge. A small dead water zone exists in the immediate' ue'ighborhood behind 
the wing-wall junction which causes poor heat transfi'r at that location. However, in 
the downstream of the wing, heat transfer is increased remarkably as compared with the 
plane chamiel flow. As such, even for a very long channel (X = 10.28), the ('nhancement 
of the heat transfer at the exit of the channel is more than 4() perc('ut. TIk^ enhancc'ment 
is not so pronounced when the effect of a punched hole beiu'ath th(i wing is takc'u into 
account. Due to the downward normal stream at the cross-plaiu'.s (as shown in Fig.4.7), 
the strength of the longitudinal vortices is reduc(Hi to a great <'xt('nt although a spiraling 
flow pattern still exists. The improvement in tlu' heat transfi'r cof'fficient is r('lativ('ly less 
than that of the case without any punched hole. How('ver, at tlu' (!xit of the chaniK'l, the 
improvement of heat transfer is about 29 percent as compared with th<^ plant' channel. 

Figure 4.9 shows the lateral distribution of local heat transfitr cot'fficituits (local 
Nusselt numbers) at different axial locations bt'hind tlu' wing on tln^ bottom plate of 
the channel. As one moves downstream, the peak of tint local Nusstdt number rt'duces 
and the Uiteral extent of the heat transfer (uihauctuiu'ut widt'us slightly. 

The effect of Reynolds number on combined spanwise avt'rage Nussc'lt nunilx'r is 
clearly evident from Fig. 4.10. Higher Reynolds numlx'r signifit's grt'att'r mass flow rat(' 
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Fi^iut. 4.7. Effec t ef stamping on cross-stn^aiu velocity vc'ctors at diffeu’ent axial locations 
in the channel with delta wing as th(^ ob.stacle 
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Figure 4.8: Effect of stamijing oil the distribution of combined spanwise average Nusselt 
number in the channel; simultaneously developing flow 
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Figure 4.9: Lateral clistributiou of local Nasst'lt uumber at differcut axial locatioiLS in 
the chamiel 
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and consequently liigln'r vortt'x strength is obtainexi ensuring be'tter heat removal. At 
the exit of a rnoeie'rate'ly long chamwl (A" = 8.25), Nu,„ for R.c = 1000 is 37 percent 
higher than that for i?e’=500. Howe^vea', at the' same' letcation, improve’rnents of 69 and 
96 percent e)ve'r the' e-ase' e)f R.(> = 500 are ol)se!rved for the Reynolds mmibers of 1500 
and 2000 respectively. In each case, the Nusselt mmiber distribution for plane channel 
flow has been she)wn by the dotted lines. In this study, tlie length of the channel is 
well below the thermal entry length. What we have illustrated here is the influence of 
Reynolds numl)er on the enhancement of heat transfer in a thermally developing flow. 
It may be mentioned that tlie fully developed Nusselt number (in our case Nu = IiH/k) 
for a channel flow with uniform wall temperature boundary condition is 3.77 . 

Augmentation of heat transfer is ol)tained at the price of increasing pressure drop. 
As we have mentioned earli('r, the increase in pressure drop is relatively less when wing- 
type vortex generators arc^ used for enhancing heat transfer coefficients. In order to 
show the aforesaid performance of wing-type vortex generators, the combined spanwise 
average friction coefficient has been defiiKHi as 




Figure 4.11 shows th(^ effect of varying tlie vortex generator’s angle of attack keeping 
the size constant. Increasing the angle of attack has the effect of increasing the vortex 
circulation wliich increases resistance and a higher value of combined spanwise average 
friction coefficient is obtained. We may mention here tiiat our friction coefficient is not 
similar to Darcy’s friction factor, rather, it is average Fanning’s friction factor. However 
at the exit of the channel, (C/ x Re) for = 25" is 60 percent more than that for the 
plane channel and (C*/ x Re) for /? = 35" is 12 percent more than that for ^ = 25". We 
have considered thermally developing and hydrodynamically developed flow herein. For 
the plane channel, {C j x R,e) value remains unchanged (=12) throughout, which is the 
exact solution for plane Poiseuille flow. This observation is a check for computational 
accuracy as well. 

Figure 4.12 shows the distribution of combined spanwise average friction coefficient 
{Cf X Re) for three different cases. The plot clearly shows the effect of stamping on 
the distribution of the (C/ x Re). Here we have considered a simultaneous^ developing 
flow. For a vei'y long channel, at the exit (AT = 10.28), the increase in (C/ x Re) for 
the case with a built-in wing-type vortex generator ( without stamping )j)ver that for a 
plane channel is 64.71 percent, whereas this value of friction coefficient {Cj x Re) is only 
49.44 percent more than the plane channel value when the effect of the punched hole 
under the wing is taken into account. The reason for tliis reduction in friction coefficient 
is attributed to the spiraling flow with relatively less vortex strength. Strength of the 
vortical motion is reduced due to the intei'ference with the downward normal stream 
created at the location of the punched hole wliich is also evident from Figure 4.7. 
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Figure 4.10: Effect of Reynolds nuinber on the distribution of coinbiiu'd sp;iuwis(' avc'rnge 
Nusselt luunber in the channel with and without the delta wing 
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Figure 4.11: Effect of angle of attack of delta wing on the distribution of combined 
spanwise-average-friction-coefficient in the channel 
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Figme 4.12. Effec t of .stairiping on th(> cii.stributiou of eoiubinc'd spauwiso avc-ragc' friction 
coefhcieiit in the cfiannel; simultauc'ously dcwciloping How 
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Figure 4.13 show.s tlu' ('fleet of Reynolds nuinlx'r on (C; x Re) while keeping thc! 
pararrK't.c'r.s r('lat('(l to_ehaunel and Mk' vorR'x g('nerat,or constant.. At thc (ixit of thc 
channel {X = M-24), {C / x Re) foz_/?:e=l()()() is about 42 pc'rccnt grcab'r than that for 
R.e=[)0{). At th(' sanu' location, (C; x R.e.) for 7?e=150() and i?,c=2()0() arc about 76 
percent and 101 percc'nt gr(^at<'r t lian that for /?e=5()() respectively. A liigher mass flow 
rate increases th(' strength of tlu^ vortices winch eventually results in a higher drop in 
skin friction. 


4.3 Performance of Delta Winglet-Pair Type 
Vortex Generators 

As we hav(' di.scussed earlier, tlu' flow in a channel with an attached delta wing shows 
no trailing edg(' vortic('s. For th(' configuration of int(!rest, the main difference between 
the wing and wiugh’t is that the wingltd. has a fre('. trailing edge whereas the wing has 
no trailing edge (wing si)au is attached to the i)late). 

Figure 4.14 shows ('ross-str('ani velocity vt'ctors at various axial locations of the 
channel with a luiilt-in wingl('t-j)air. For the winglet-pair the vortices do not have to 
turn as abruirtly near the? trailing edge as they do for the wing. However, in order to 
have a lucid understanding of the structure', secondary flow, various magnified scales for 
the average axial vcdocity hav(' b('.(,'n used for different cross-sti'eam planes. 

Normalized strc^ainwisc' vorticity contours for the above mentioned case have been 
shown in Fig. 4.15. Tlu' dottcxl liiu's indicat(^ negative vorticity. Peak vorticity occurs 
in the vortex centers and m^ar th(' symmetry sides (side walls) where the symmetry 
boundary conditions lead to rapid change in the direction of flow of the secondary 
vortices. 

Figuz'e 4.16 shows the isotherms at different cross-planes in a channel in presence of a 
delta winglet-pair. The vortex generators influence the temperature fleld strongly. The 
spreading of isotherms over a wider region of the cross-section implies that the bulk 
temperature is higher in thc^ downstre^am from the trailing edges of the winglet-pair. 
The importancfj of the secondary flow and the liigher heat transfer is evident from this 
picture. 

Figure 4.17 shows the influence of angle of attack of the winglet-pair on combined 
spanwise average Nusselt number in the channel (equation 4.1). At a nondimensional 
distance of 4 from the inlet, for an angle of attack of 20^' , we observe an enhancement 
of 30 percent in the combiiwai spanwise average Nusselt number over the case of a plane 
channel. Now, at the same location, for an angle of attack of 26" , an improvement of 
about 6 percent in Nu^a over the; case of the 20" angle of attack is obtained. However, 
at the same location, an improvement of 12 percent over the case of 20" angle of attack 
is discerned for 0 — 32". The; wiuglets with higher angle of attack produce vortices with 
higher strength, which in turn results in improved heat transfer. 
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Figure 4.16; Isoth('.rrn.S' at. difftu'eiit. eross-planes in a channel in presence of a built-in 
winglet-pair 
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Figuro 4.17: Effect of angle of attack of wiuglet-pair ou th(^ distributioa of conibiiu'cl 
spaiiwisfc-avt,iage-Nu.sselt-uumber iu the channel (level()[)('<l profile at the inlet 



4.4. COMPAIUSON WITH EXPEIUMENTS 


81 


Figure 4.18 shows Mu' ('fh'ct, of wiuglef-pair on disfrilmtiou of in a channel 
for various Reynolds nuinhers, Wingh't.s with high('r Reynolds number induce higher 
vortex stnuig t h a nd t h<'r('l)y iiigher iK'at, transh'r is brought, about. For all the Reynolds 
uumbtus, tlie Nu,„ distriluition has Ixuai compared with the corresponding plane channel 
value. At a nondiiiK'usional distance' of 2.5 from the inlet, for a Reynolds number of 
2000, we obs('rv(' an (uihance'nK'nt of 30 ix'rceuit in the' Nus„ over the corresponding 
value for a chamK'l without any obstacle'. 

Figure 4.19 shows the e'ffe'ct of varying die angle of attack of the winglet on combined 
spanwise ave'rage' skin friction coe'fficient (C/ x 7?.e;), givem by eepiatioii 4.2, wliilc keeping 
the size constant. Incre'asing the' angle' e)f attack has the effect of increasing vortex 
strength, which in turn ine-re'ase's re'sistance' anel conse'eniently a higher value of combined 
spanwise' ave'rage' friction e-e)e'fiici('nt. is e)btaine'el. At the' exit eif the channel, {Cj x Re) 
for 3 = 32" is ne'arly 8 jie'rce'nf more' t.han that feir ,3 = 20" and (C/ x Re) for /3 = 26" 
is about 3 pe'i e-e'nt highe'r tlian that for ji = 20". 

Effect e)f Re'yne)lds numbe'r on the* distributiem eif (C; x Rxi) for the built-in winglet- 
pair has lie'on shown in Fig. 4.20, Ine-re'asing (C; x Re) feir higher Reynolds number 
is a e:e)nseeiue'ne-e' e)f iimre'a.se'd vort.e'x st.re'iigth for liiglmr fleiw velocities. For Reynolds 
number of 2000, the' {C/ x Re) is 82 iie're-e'ut me)re than that of Reynolds number of 500 
at the exit e)f the' e’hanne'l. 


4.4 Comparison with Experiments 


The moded valieiatiein is iie'rfbrme'd tlireiugfi comparison with some published experi- 
mental results. Le)e-al Nussedt numlieirs aleing the centerline of the bottom plate with a 
delta-wing at an angle of attack of 20" were calculated. Reynolds and Prandtl numbers 
for this computatie)!! are. 1815 anei 0.7 respectively. Here local Nusselt number values 
are evaluated on tlie basis of entry temperature of the incoming stream as 


"k[(g),.o,..o].(H/k) 


(4.3) 


Figiu’e 4.21 shows that the computed values of local Nusselt numbers (with 
weighted averagt^ scIkuxk^) comparx; favorably with the experimental results of Fiebig 
et al. (1986). In the ('xix'rinumt (Fiebig et al, 1986), the channel walls had punched 
holes on them and tht' (^xix'riiiK'utal results are closer to the computed results for the 
case with stamping. Tln^ comixutational rtisults due to QUICK discretization scheme 
(Leonard, 1979), tiavci also agnsc'd (dosely with the experiment with the exception near 
the exit of tlie channel. Tfw'. lo(;al Ik '. at transfer coefficients were determined by unsteady 
liquid crystal tliermography Measurement error is indeed minimal in this process. The 
RSS uncertainty dx'seribed by Moffat,(1987) was below 4.3% for local heat transfer coeffi- 
cients. Howev('r, t.h<' small discrepancy Ixdwcxm tfie experimental and numerical results 
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Figure 4.18; Effect of wiuglet-pair ou combiiied spauwise av<'rag(; Nuss<'lt iiumfx'r dis- 
tribution of tfie chaimel for different Reynolds imrnl)crs 
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could bo at tributed to ])bvsioal boundary ooudit.iou on th(^ channol walls. Duo. to very 
small thiokiH'KS (H.4 x 10 ' ni) of t.lK' fin-plat c' (us('d iii o.xporimont), tho axial heat cou- 
duction was not. goo<l and it was not possible' to maintain porfoct isothermal condition 

on the chaniK'l walls. 

Overall boat, transh'r i)orformano<' data is validated tlirough comparison with respect 
to mean-Coll)urn-faotor (j) using d('lta-wing for three different angles of attack. The 

eriinental n^sults wore obtained from Fiebig el al (1991). In the experiment, the 
Serence plate-fin area started from tlu' location of the base of the wing; its longitu- 
dinal extent was 7.5 wing chord lengths and its lateral extent was 4 times the wing 
span Identifying this /,on<' on tlu* bot tom platt' as the, reference area, experiments were 
conductc'd in a idaiu' chaum'l and in a channol with a built-in delta wing at different 
angles of attack. Figniv 4.22 shows t lu' nu'an-Colburn-factor(j) versus Reynolds number 
for a delta-wing of asix'ct rat io 1.25 at various angles of attack. The numerical model 
and' the experiment corroborate' wit h ('ach otlu'r reasonably well. The curves can be 


represente'd l>y i~( , 


wIk'TC' u uihI (\i Hr(' (‘oiistauts. 





Figure 4.21: Compm-isoii of computed results with expmiuKmtal olisvrvatiou 
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Figure 4,22: (Comparison l)(‘tw(‘(ai (■omi)u(.at.i()nal and expcriiiieiital results for Colburn 
factor as a fmuiJon of H<*ynol(ls numlxT ai diflenajit angles of attack 



Chapter 5 


A Second Law Based Optimization 
for the Choice of Vortex Generators 
in the Laminar Flow Regime in the 
Laminar Flow Regime 

"At Jir.st ulaiicv eiigineeiing thermodynamics may 
appear to be a dead discipline, but closer examina- 
tion reveals a remarkably lively field of study. Its 
(nrionis have recently been reformulated in terms of 
heat transfer in.Htead of mechanics. New graphic 
const ruet ions are being used to summarize proper- 
ties. And exergy analysis and refrigeration tech- 
niques are active pursuits for many engineers" 

A. Bejan, Mechanical Engineering, p58, May 

19HS. 


5.1 Introduction 

Cliaptor-4 gives a (jaant.it.at iv<' ix'rfonuaucc data for a delta wing and winglet-pair 
in a chauu(4 wi(,h (liflerent angh'S of attack and a wide range of operating conditions. 
Howevcir, the ixuddnuaiicf- of wings and winglets can be compared by taking into account 
two i)aram('t<'rs: the augiueut.at ion of Iwuit transfer and the associate flow losses. In the 
case of conv<u:t,iv(‘ heat, t.ran.sf(*r t.wo types of losses of useful work is encountered; one 
due to fluid friction and ot her <lu<“ t.o heat transfer across finite temperature gradient. 
Both th('.se pheiioiiKaia are inanih'stations of irreversibility and therefore for a complete 
evaluation of h<>at, t.ransfer augmentation tedmique a detailed thermodynamic analysis 
is requin'd. Thus, t he i)r('S('nt chapt.r'i' incorporates aii. detailed analysis on the influence 
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of vortex generators (wings/winglets) on entropy generation and irreversibility. An 
evaluation of the exergy exchange process from the standpoint is called second law 
analysis (Bejan, 1977; Bejan,1979). Based on an evaluation of irreversibility, this study 
is able to draw some conclusions about the efficient use of thermal energy. 


5.2 Irreversibility Approach to Heat Transfer 
Process - A Brief Review 

The concept of irreversibility can be employed as a general criterion for designing 
various thermal systems in order to minimize the wastage of usable energy (McClintock, 
1951; Bejan, 1977; Bejan, 1979; Nag and Mukherjee, 1987). Irreversibility analysis 
through the use of Second Law of Thermodynamics considers the utilization of available 
energy (exergy) in contrast to the usual design procedure aimed at maximum transfer 
of heat in a heat exchange process. 

McClintock (1951) is one of the pioneers who have used the Second Law of Thermo- 
dynamics in heat exchanger analysis. He has studied the production of irreversibility 
in a heat exchanger designed to exchange a specified amount of heat between the flow 
streams. 

Irreversibility analyses by Bejan and Smith (1975) have successfully optimized the 
heat exchanger passage area along the flow path of heliiun vapor. 

Bejan (1977) has introduced the concept of designing heat exchangers for fixed or for 
specified irreversibility in contrast with traditional design procedure based on specified 
amount of heat transferred. This approach is particularly useful in design of cryogenic 
systems where conservation of available refrigeration is of prime consideration. 

Irreversibility analysis for fundamental convective heat transfer has been carried out 
by Bejan (1979a). Four fundamental flow configurations have been illustrated: forced 
convection in a round tube, flow over a flat plate, single cylinder in a crossflow and 
flow in the entrance region of a rectangular duct. The expressions obtained from this 
analysis serve as the basis of analyses for heat exchangers. 

For insulation systems where the ratio of absolute temperatures across the insulation 
is liigh and where the conductance of the insulation is already reduced to the minimum 
acceptable economically, Bejan (1979b) has prescribed the temperature profile across 
the insulation that conserves the available refrigeration. 

Sarangi and Chowdhury (1982) have analyzed counter flow heat exchangers to ac- 
count for the entropy generated due to axial conduction and have tlerived an expression 
for optimum wall conductivity. 

Nag cind Mukherjee (1987) have shown that for a heat transfer process tlu're (exists 
an optimum ratio of heat transfer to [dumping power. Simply maximizing tliis ratio is 
not often a good decision, since in that cas<'. tlu; (uitropy generated may be far from 
the rninimiun possible and a larg<! amount of th(^ available (uiergy inny thus Ix^ irre- 
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The first term on the right-hand side of (5.2) may be written as . 


^(g.VT) = ^ 


dx J ^ [dy J [dzj 


The second term on the right-hand side of equation (5.2) is expanded as 


(5.3) 


2 ^ f dui 

{cT-.VV) = (p+-axA)%( — 


( dui duj\ dui 
^ i dx 


dxj 


or, (a-.W) = (p-h-/iA)A-/i 


+ 


+ 


+ 


' du dw\ du 

^dz dx j dz'^ 
' dv dw\ dv 
, dz dy j dz'^ 


du du [ du dv \ du 
^dxdx \dy dxj dy 
' dv Su'i dv dv dv 
dx dy j dx dy dy 


^ dw du\ dw 
, dx ^ dzj dx 


'dw dw dw dw 

dy ^ dzj dy dz dz 


After invoking the incompressibility condition, equation (5.4) can be written as 


{a : Vy) 


E 

T 



+ 


' du dv^ 
dy ^ dx, 


+ 


' dv dw' 


( du dw\ 
^ \ dz dx ) 


(5.4) 


(5.5) 


Use of equations (5.3) and (5.5) in equation (5.2) and nondimensionalization yields ; 


Ns. = 


S.H- 


(r - 1)'^ 


+ 


+ 


[0(r-l) + ip- 

Ec.Pr 

[%-l) + l]- 

fdu dV' " 



+ 


'd9' 

dY. 




dx \dY, 


+ 


'dw' 

dZ 


dV dW' 
dY ^ dx) ^[dZ^ dY 


(dU dW' 
'^\dZ^ dX 


(5.6) 


As expected, the irrev('rsibility indicator AT.s'rj contains two additives parts, on<! due to 
conduction in the presence of non- zero temperatiu'e gradicuit, and the other accounting 
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for viscous dissii)atiou of iiu'chanical ])owor throughout, the flow. However, if the' velocity 
and te'inperature fields an' comijlete'ly known, we can find out tlie value's of ne)neiimen- 
sie)nal entrerpy geuieu'ation jx'r unit volume (Ns^) at each cell e>f the fle)w dennain. After 
evaluating this ejuantity, we' inte-grate' TV.s's over the entire volume' to get nondimensional 
vedumetric entropy generation in the cliannel as ; 

Ns = IJJ Ns 3 dXdYdZ (5.7) 

5.3.3 Merit Function 

From this analysis, it is possible to evaluate the rate of energy transferred usefully as 
well as destructiem of exergy due to irreversibilities^ _ 

If Q is the total rate of heat transfer, then Q =hc{T^. - Tb) 2BL 

or, Q = 2{B/H) N^c kL (r„, - Tf.) (5.8) 

The equation (5.8) can be rearranged as : 

Q = 2{B/H) kimT, (5.9) 

The rate of exergy transfer accompanying energy transfer at a rate of Q is given by 
Moran (1982 ) as 

0. = q[i-|]=q[i-;] (5-10) 

where T^ , the exergy reference environment temperature has been considered as the 
ambient temperature Too ■ The wall temperature Ty^ has been considered as a suitable 
temperature at the surface where heat transfer takes place. If S is the total rate of 
entropy generation (dimensional), the destruction of exergy (Gouy-Stodola theoiem) 
is 

I = TJ^- Ty,S (5-11) 

r 

As such, the total rate of entropy generation, 5, may be written as 

J d{HX) d{HY) d{HZ) (5.12) 

Invoking equation (5.7) in equation (5.12), we obtain 

S = HkNs (5-13) 

Finally, equation (5.11) and equation (5.13) will yield 

/ = - Ty, Hk Ns 
r 


(5.14) 
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A merit function is defined as the ratio of exergy transferred to the sum of exergy 
transferred and exergy destroyed 


M = 


Qa 

<3a + / 


(5.15) 


Substituting for Qa 
obtain, 


from equations (5.9) and (5.10), and "I" from equation (5.14), we 

, , 2(B/ H) (r — 1)" Nuc (&w ~ ^b) ('5 16') 

“ 2(B/H) (r - 1)2 Ma (0w -A) + f Ns 


This merit function is now evaluated for various flow parameters in a channel using 
delta wing or winglet-pairs as vortex generators. However, it might be observed that 
equation (5.16) is another form of second-law efficiency (Moran, 1982 ). The merit of 
the merit function lies in its simultaneous accountability of exergy transferred and its 
destruction caused by irreversibilities associated with exergy transport and momentum 
transport. Irreversibilities due to external interaction and internal dissipative effects 
are together taken care of by this parameter. 


5.4 Results and Discussions 

This chapter contributes a quantitative comparison between the performance of a delta 
wing and that of winglet-pair on enhancement of heat transfer. It will be of interest to 
accomplish the comparison by taking into account the irreversibilities associated with 
the heat transfer process. 

Figure 5.1 illustrates the relative heat transfer performance of delta wing and winglet- 
pair configurations with other geometrical and flow parameters unchanged. The main 
difference between the wing and winglet is that the winglet has a free trailing edge, 
whereas the wing has no trailing edge ( wing-span is attached to the plate ). However, 
delta- wing produces streamwise vortices with higher strength, which brings about better 
improvement in heat transfer as compared with built-in winglet pair. In a moderately 
long channel (X = 8.4), the enhancement in combined spanwise Nusselt number for 
a built-in delta-wing (/3 = 26°, A=1 and Br=0.42) at the exit of the channel, is more 
than 34 percent than that of a plane channel. This enhancement for the winglet-pair 
at the same location is about 14 percent. Therefore from heat transfer point of view, 
delta- wing is found to be more effective than the winglet-pair. 

_ Figure 5.2 shows the distribution of combined spanwise avertige friction coefficient 
(Cf X Re) in the channel for the wing and winglet-pair and compares them with the 
distribution of (Cy x Re) in a plane channel. As mentioned earlier, delta- wing gcnier- 
ates streamwise vortices with higher strengths. Increased frictional losses are due to 
steeper velocity gradients at the walls. These losses rise witli vortex strcuigth, and as a 
consecpience, (C/ x Re) is higher for vortices with higher strength. At the chaniK'l exit 
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Figure 5.2: EflFect of type of obstacle oa the distributioa of combined spanwise average 
friction coefficient in the channel; developed profile at the inlet 

[X = 8.4), (C/ X Re) value for a delta- wing is about 79 percent more than for that of 
the plane channel flow. For winglet-pair, the increase in (C/ x Re) value at the same 
location, over the case of a plane channel flow, is nearly 65 percent. 

Figure 5.3 shows the total entropy generation in the channel at various Reynolds 
numbers (in the range of 500 — 3000) for two different cases, namely, with a built- 
in delta wing and winglet-pair. It is evident that for all the Reynolds numbers, the 
entropy generation in the case of a built-in delta wing is much more than that of a delta 
winglet-pair. With increasing Reynolds number, the delta-wing generates vortices of 
higher strengths. This culminates in steeper velocity gradients at the wail. As a result, 
frictional losses become high and consequently, the total generation of entropy becomes 
much more pronounced. 

Figure 5.4 shows the variation of merit function, M with Reynolds number (in the 
range of 500 - 3000) for the built-in delta wing and winglet-pair. The delta winglet-pair 
shows better performance than delta wing with regard to merit function. For the case 
of heat transfer in a channel with a built-in winglet-pair, the incn-ase in irreversibility 
associated with the increase iii exergy transha’ occurs at a relatively lower rate as com- 
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re 5.4: Effect of type of obstacle on the variation of merit function with Reynolds 
ber 

i to the case of a built-in delta wing. Although with regard to augmentation of 
transfer in the channel, the built-in delta wing shows better performance, it can be 
linly argued that the built-in winglet-pair is more effective with respect to efficient 
)f energy. 


Appraising Remarks 

/ith above mentioned events in mind, we are now in a position to consider the 
tical application of this augmentation technique. The flow configuration resembles 
ly a single element of a cascade in fin-tube or plate-fin crossflow heat exchangers 
a row of vortex generators. In practice, fjr such an exchanger tht; channel may 
be allowed to be so long (X > 8.0), i.e., another row of vorttix generators may 
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not bo alkwod to b<' so long (A' > 8.0), i.o., another row of vortex generators may 
be pnnch('(i out or inountc'd bcdore tliis k'ligtli is n'ached in the downstn'arn direction. 
In that event, at about A" = 4, (refer to Fig. 5.1), augmentation of heat transfer due 
to winglet-pair is also quite substantial. Moreover in the case of winglcts, separation 
bubbles are not formc'd near the wing-plate junction. It has also been seen that the flow 
loss (corresponding to the combined spanwise average skin friction coefficient) due to the 
winglet-pair is less than that due to the wing. Besides, due to formation of trailing edge 
vortices for the free trailing edge of winglet, the small zone of poor heat transfer as it is 
observed with the wing can l)e avoided. Thus, based on tliis and a comparative study of 
second-law efficiency, the use of winglets appears to be a more attractive augmentation 
technique'. 



Chapter 6 

Turbulent Flow and Heat Transfer 


6.1 Introduction 

It may be mentioned that the study of laminar flow in the present work is not for 
computational simplification. Usually, the fin spacing is so small and the mean velocity 
range is such that the flows are often laminar (Kakac et al, 1981). However, for some 
special applications, the velocity becomes such that one may encounter turbulent flows 
in plate-fln type heat exchangers. Therefore, a detailed analysis of turbulent flow in 
an element of such heat exchangers is necessary. It may be mentioned that in plate- 
fln heat exchangers no holes should exist underneath the vortex generators in order to 
avoid mixing of the hot and cold streams. It is conjectured that in a laminar flow the 
coherent motions are dominant over the chaotic behavior of the flow. In the transitional 
flow chaotic and coherent behavior are equally dominant. When the chaotic behavior 
becomes dominant in relation to the coherent behavior, the flow is turbulent. 

According to this statement the turbulent motion is chaotic and irregular. However, 
if it is entirely irregular or chaotic, turbulent motion will be inaccessible to any type 
of mathematical treatment. Instead, it can be assumed that the irregularity associated 
with turbulence is such that turbulent motion can be considered as an irregular condi- 
tion of flow in which various quantities, viz., velocity components and pressure show a 
random variation with time and space so that the statistical average of those quantities 
can be described. 

The turbulent flow is dissipative in nature. Therefore if there is no external source 
of energy for the continuous generation of turbulence, the turbulent or irregular motion 
within a flow will decay. It is postulated that the fluctuations inherently come from 
disturbances and may either be viscosity damped or may grow by drawing energy from 
the free stream. At a Reynolds number less than critical, the kinetic energy of flow 
is not enough to sustain the random fluctuations in the face of viscous damping and 
in such cases laminar flow continues to exist. At somewhat higher Reynolds number 
than the critical Reynolds number, the kinetic energy of flow supports the growth of 
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fluctuations and transition to turbulence is induced. 

However, the self-sustaining character of turbulence in a simple shear is that the 
turbulent velocity fluctuations in the direction of mean velocity gradient promote a 
growth in shear stress which in turn, serves to augment the intensity of streamwise 
fluctuations. Pressure interactions deflect some of this energy to fluctuations in the 
direction of the velocity gradient and thus the sequence repeats itself. Our aim is to 
mimi c these processes in terms of time-mean properties. As such the most popular 
approach to resolve the issues of turbulent flow is to solve time-averaged (or Reynolds- 
averaged) equations and model the closure equations through some hypothesis. This 
approach has been widely used in engineering, but suffers from being problem dependent 
and non-universal. As we have seen, an extensive literatiue exists on the performance of 
the various turbulence models. The most popular turbulence model is the two-equation 
model. At the practical level, satisfactory results are obtained by making use of the two 
equation turbulence model due to Launder and Spalding (1974). 

This chapter presents numerical simulation of flow and heat transfer in the configu- 
ration of interest (Fig. 3.1) for the turbulent regime. The governing equations and the 
boundary conditions have been discussed in section 3.2.3 and 3.2.4 respectively. Now 
we would like to highlight the basis of derivation of turbulent kinematic viscosity. 

By analogy to molecular viscosity, we can write 

r'i = VsL (6.1) 

where, Vg is the turbulent velocity scale and Ig is the turbulent length scale. 

On the other hand, k, the turbulent kinetic energy is given by 

k = + ^ (6.2) 

Therefore a meaningful scale for velocity fluctuation is k^^-. The length scale is related 
to the rate of dissipation of tiubulent kinetic energy, e, given by 

ls = (6.3) 

€ 

Finally, the Prandtl-Kolmogorov relation is (obtained from dimensional argument as 

Ut = Cf,fc-/e (6.4) 

where is the constant of proportionality cind is equal to 0.09. The nondimensional 
form of the equation (6.4) has been used as equation (3.13) in section 3.2.3. 

6.2 Basis of the Wall-Function Treatment in Turbu- 
lence Modeling 

The no-slip condition at a rigid boundary implicates that couv('ctive transport as- 
sociated with the motion paralkd to the solid wall vaiiishes there*. Thus, very neuir 
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to the wall, the system of diffcu'cntial equations describing the mean and fluctuating 
fields become simply a function of normal coordinate, y. Therefore, with the requisite 
boundary conditions provided, the equations could be solved to obtain the distributions 
of the various dependent variables as a function of y and other relevant flow parame- 
ters. Any of the turbulence models discussed in the section 2.4.2 could serve a basis of 
generating such near wall solutions and the results obtained would be as accurate as a 
complete low-Reynolds-number near-wall treatment provided the neglect of streamwise 
variations is accomplished within the thin near- wall region {Y'^ < 11.63). This feature 
permits us to patch a tliree dimensional solution of the outer fully turbulent region to 
a one-dimensional near-wall solution in which all the near-wall influences are restored. 
The inner solution is conveniently expressed as a system of algebraic formulae, such 
as, equations (3.23) and (3.24) relating the dependent variables {U, W, k, e etc.) in 
their dimensionless form to the independent variable y (also expressed dimensioiilessly) 
and various parameters of flow. These formulae are commonly referred to as "Wall 
Functions". 

The section 3.2.4 provides with the formulae used for two layer wall function (see 
Amano, 1984) treatment. It may be mentioned that in all the formulae, the friction 
velocity, Ur is one of the most important characteristic quantity. However, in this 
section, the theoretical basis of expressing Uj in terms of other quantities related to 
turbulent flows will be discussed. 

The local equilibrium between the production and the dissipation rate of turbulent 
kinetic energy near the wall gives rise to 

pvlv' ^ = />€ (6-5) 

dy 

Just outside the viscous sublayer, the shear stress (still equal to r,,) will be produced 
entirely by the turbulent eddies (Bradshaw, 1971) so that 

pu'v' = T„, (6-^) 


From equation (6.6) we can also write 

, du 

u'v' = T^,lp = Vt— 


(6.7) 


Invoking i/j = C^k'^/e in equation (6.7) we can write 

C^k"^ du (6 

u'v = '' 

e ay 

Substituting the value of {du/dy) from equation (6.5) in the above equation, we get 
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or 


(-W)' = C^k- 


or 


Invoking Ur = 



p 


In nondimensional form, equation (6.9) may be written as 




(6.9) 

( 6 . 10 ) 


The friction velocity is of paramount importance in deriving out the nondimensional 
form of the governing equations and the boundary conditions for turbulent flows. 


6.3 Present Investigation 

An effort has been made to investigate numerically the turbulent flows and associ- 
ated temperature fields in channels with longitudinal vortex generators. The present 
investigation serves the following purposes: 

(a) validation of our model through comparison with the experiment results of 
Pauley and Eaton (1988 a,b) and Eaton (1993). 

(b) prediction of the heat transfer and flow losses in the configuration of interest. 
We envisage to predict the distribution of spanwise-average-Nusselt number 
and the distribution of spanwise-average-skin friction coefficient on the bot- 
tom plate of a channel with built-in vortex generators. It is implicit that the 
flow regime is turbulent. 

6.4 Comparison of Numerical and Experimental Re- 
sults 

6.4.1 The Test Section and the Computational Domain 

Experimental data for three-dimensional tiu'bulent flows in a rectangular channel 
containing longitudinal vortex pairs have been presented in Pauley and Eaton (1988 
a,b) and Eaton (1993). Figure 6.1 shows the schematic of the experimental facility. 
This facility has a 200 cm long test .section with a 13 cm x 61 cm cross-section and is 
operated at a normal free-stream velocity of = 16 m/s. A pair of delta winglets is 
used as vortex generators and is located at 53 cm dowustn^am from the inlet. Tln^ fUigk's 
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Figure 6.1: Schematic of the test section of Pauley and Eaton (1988 a,b) 

of attack of the' vortex generators can be varied. The case with the winglet height of 2 
cm and a length of 5 cm with 18" angle, of attack, which has been referred to in Pauley 
and Eaton (1988 a,l>) and Eaton (1993), is chosen for the present simulation. Pauley and 
Eaton (1988 a,b) and Eaton (1993) have conducted measurement at four cross-sections, 
namely, x=66, 97, 142 and 188 cm (see Fig. 6.1). The section x=97 cm is selected as 
inlet of the computational domain used for the comparison. No experimental data is 
available for the turbulent variables at the axial station of .'C=66 cm. The computational 
domain ends at the <!xit of the channel (i.e., at x=200 cm). The measured data of Pauley 
and Eaton (1988a) for u, v, w, k and e at the cross-section a:=97 cm are deployed as 
inlet profiles of these variables for the computational domain. 


6.4.2 Comparison and Discussion 

The flow is simu-lated for a Reynolds number of i?eH/2=67000, based on half height 
of the channel. Air has been used as a working fluid; hence Prandtl number of all 
computations is 0.7. 

The comparison of secondary velocity vectors at three different planes in the channel 
can be observed through Fig. 6.2. A reasonably good agreement is confirmed between 
the experimental data and computed results. Both in the experiment and in the com- 
putation, the vortex transport and decay in the strength of secondary ow in various 
cross-sections are discerned. 

Figure 6.3(a) shows isolines of the streamwise velocity at three different cross-stream 
planes of the channel. The boundary layer is thinned in the downwash region and seems 
remain unaffected in the upwash region. The vortices diverge as they move downstream. 
All such salient features are captured in computational results w c are i us ra e in 

Fig. 6.3(b). 
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Figure 6.4(a) shows the experimental values of isolines for the turbulent kinetic 
energy at two different cross-sections in the channel. In the downwash region, the 
boundary layer is thinned by the main vortex and the location for the peak value of 
turbulent kinetic energy comes closer to the wall. In the upwash region, the iso-kinetic 
energy lines are more uniformly placed. Although the k-e model formally employs some 
restrictive assumptions, the results obtained from the present computation (Fig. 6.4(b)) 
seem to match surprisingly with the experiments. 

An engineering computation of turbulent flow concerns determination of mean axial 
velocity profiles in a flow field and the effect of time- dependent fluctuations on them. 
Figures 6. 5 and 6.6 compare the experimental and computed axial velocity profiles 
(normalized) at different axial distances for five spanwise locations. In Fig. 6.5, the five 
spanwise locations at z=2.4, 4.3, 6.2, 8.1 and 10.0 cm represent the central core region 
(A), downwash region (B), region of vortex core (C), upwash region (D) and outer vortex 
region (E) of the cross-section at x=lA2 cm while in Fig. 6.6 the same spanwise locations 
stand for the similar five regions of the cross-section at x=188 cm. In the central region 
(A) of the cross-section at x=142 cm where the boundary layer is tliinned significantly 
by the influence of strong downflow (refer to Fig. 6.2), the velocity profile shows a 
usual two dimensional distribution. However, in this region the agreement between 
the computational and experimental results is indeed good. In the upwash region (D) 
the streamwise velocity profile is not so correctly reproduced by the simulation. The 
departure of prediction from the experimental result goes up to 6 percent. In other 
regions of the cross-section at x=142 cm, a reasonably good agreement between the 
experimental results and the prediction is observed. Neiirly similar trend of comparison 
between the experimental observations and the predictions is illustrated in the plots at 
x=188 cm. 

Figure 6.7 compares the computed and measured values of tiubulent kinetic energy 
distribution in the y-direction for different z locations at an axial station of x=188 cm. 
In some cases, (at z=6.2 cm), the computed results marginally overpredict the results of 
Eaton (1993). By and large, the comparison implicates acceptable agreement between 
the prediction and the computed results. 

6.5 Prediction of Pertinent Performance 
Parameters 

The predicted static pressure distribution at different axial locations in the channel 
for a Reynolds number of Re — 5000 is shown in Fig. 6.8. It is noteworthy that the 
longitudinal vortices persist in the downstream section influencing th(^ e.ntire velocity 
field. The basic free-vortex like structure is also observed in case of th(^ turbulent 
flow. The location of the minimum pressure is at the vort(!x ccuiter. The minimum 
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Figure 6.4: Isolines for turbulent kinetic energy at different cross-sections for 
Ren/2 = 67000 








U/Uo 


ReH/2 = 67000, /3 = 18 


Figure 6.5: Stremnwise velocity profiles at axial station ;/:=:142 cm for Rci 
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pressure zone moves downwards witli increasing distance behind tlie trailing edge of the 
winglets. Figm'c 6.9 confirms a similar trend of static pressure for a Reynolds number of 
Re = 15000. Figure 6.10 shows the isolines for turbulent kinetic energy at two different 
axial locations in the channel for a Reynolds number of Re = 5000). The production 
of turbulent kinetic energy is maximum in the near-wall- region where the shear stress 
is also maximum. However, due to the influence of strong secondary flow another high 
kinetic energy zone is observed near the vortex center. 

For the turbulent flows the transport rate for momentum is high and the vortex 
generators are rather small in size with a view not to enhance much of pressure penalty. 
Hence the height (b/2) of the winglets for such computations are very small compai-ed 
to the channel height ((b/2)/H=0.3). 

The heat transfer and skin friction coefficients are the two most important perfor- 
mance parameters concerning any heat transfer process. In the following text, we shall 
discuss about these two parameters. 

In order to predict the heat transfer performance at the wall, the local Nusselt 
number may be defined as 

hH 


N'U,(^x;i) 


k 


or. 


or. 


or. 


Nu{x, z) = h{Tu, — Tb) 


Nu{x,z) = 


u 1 


1 


UoH 


ot pcp Uo{T.^, - Tb) V 
q^,PrRe 


pCp{Tx. - Tb)Uo 


Nu{x, z) = 


Qw.hP'^ R>^ 


(6.11) 


where. 


0^ 

~ pCy(T„,- roo)17o 

Finally, the spanwise average Nusselt number on the bottom fin surface may be ex- 
pressed as 

, = 4 [ ^ Nu{x,z) dz (6.12) 

_ JL> Jo 

Figure 6.11 shows the Nusi distribution on the bottom plate for the Reynolds number 
of 5000 and 15000. In both the cases, the turbulent intensity is 10 percent at the inlet 
of the channel. Significant enhancement is observed for the cases with built-in delta 
winglets over the corresponding cases without any obstacle. For example, in the case 
of Reynolds number of 5000, the enhancement is about 16 percent at the exit of the 
channel. It may be worthwhile to mention here that the size of the winglets used in 
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Figure 6.8: Isoliues of static pressure at different cross-sections for /2e=r)()G0; trailing 
edge of the winglet is located at X=2.5 from the inlet 
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Figure 6.9: Isolines of static pressure at diflPerent cross-sections for i2e— 15000; trailing 
edge of the winglet is located at X=2.5 from the inlet 
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Figure 6.10; Isolines of turbulent kinetic energy at different cross-s(;ctions for i?e=r)000; 
trailing edge of the winglet is located at X=2.416 from the inlet 
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Figure 6.11: Distribution of spanwise average Nusselt number on the bottom plate for 
different Reynolds numbers 


the study of turbulent flow is indeed smaller as compared to the size of the winglets 
used for the study in laminar flow. So, the enhancement in transport rate is not as high 
as laminar flows. The peak values are brought about just inboard of the vortex center 
where the lateral divergence is most pronounced (for both the cases, it is at X=2.5 from 
the inlet). 

Figure 6.12 shows that at the location of the peak, the boundary layer is thinned by 
the strong downflow and lateral outflow of boundary layer fluid. In the downstream the 
secondary velocities decay and vortices spread apart causing a growth of the boundary 
layer. The spanwise distribution of local Nusselt number at X=4.325 on the bottom 
plate for two different Reynolds numbers are shown in Fig. 6.13. The level of heat 
transfer augmentation in the space between the vortices is high in both the cases. 
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FijjUre 6.12. Cross-stream velocity vectors at X=2.5 from the inlet of the channel for 
two different Reynolds numbers 
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Figure 6.14: Distribution of spanwise average Nusselt number on the bottom plate for 
different angles of attack 

Figure 6.14 shows the influence of angle of attack of the winglet-pair on the Nugi over 
the bottom plate of the channel for a Reynolds number of 100000. Near the location of 
the trailing edge, for an angle of attack of 12^', we observe an enhanced Nu,i (absolute 
value of 4) over the case of a channel without any obstacle. The winglets with an 
angle of attack of 45" produce vortices with higher strength, wliich in turn result in 
further improvement of heat transfer. It may be mentioned that for the turbulent flow 
the transport rate is even otherwise high as compared to laminar flow. It is evident 
in Fig. 6.14 that even for a plane channel, the Nu^i value is reasonably liigh for the 
turbulent flow. Because of this retison, for the ca.se of a built-in winglet-pair with an 
angle of attack of 45", although the enhancement in ATu,! near the trailing edge is more 
than an absolute value of 12, in terms of percentage this value is not so high (6 i>ercent). 

The local skin friction coefficient at the bottom wall may be expressed as 
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or, 

C/(x,2:) = (6.13) 

The spanwisc average skin friction coefficient on the bottom wall is expressed as 

2 

The distribution of spanwise average skin friction coefficient (C/.i x Be) for two 
different Reynolds numbers is shown in Fig. 6.15. The increase in peak values of the 
average skin friction coefficients ai-e 33 and 27 percent for the Reynolds numbers of 5000 
and 15000 respectively. 

Effect of angle of attack on the distribution of spanwise average skin friction coef- 
ficient on the bottom plate of a channel for a Reynolds number of 100000 is shown in 
Fig. 6.16. At the location of the trailing edge, the vortex strength is boosted up by the 
trailing edge vortices. The near wake regions are highly distorted by the vortices but 
in the downstream, the vortices gradually lose their strength due to turbulent diffusion. 
However, for an angle of attack of 45^', at the location of the trailing edge, the increase 
in the spanwise average skin friction coefficient (C/,i x Re) on the bottom wall is 8 
percent more than that of the case of a plane channel. To summarize, it can be said the 
transport rate is usually very liigh for the turbulent channel flows even if there is no 
obstacle present in the channel. In presence of the vortex generators, the skin friction 
coefficient (C/,i x Re)) increases by an absolute value of more than 19 but when this 
value is expressed in terms of percentage, the increase in flow loss appears to be some- 
what less. It is also to be mentioned that the size of the vortex generators for turbulent 
flows is relatively small and the transport enhancement due to the vortex generators for 
turbulent flows is not as high as laminar flows. 
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Figure 6.15: Distribution of spaiiwise average skin friction on tlie i)ottoin plate for 
different Reynolds numbers 
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Figure 6.16: Distribution of spanwise average 
different angles of attack 


skin friction on the bottom plate for 




Chapter 7 

Concluding Remarks 


The performance of heat exchanger surfaces can be significantly improved by mounting 
protrusions on tln^ surfaces. Some typical examples of surface geometries which are 
popular in different industrial applications are wavy fins, offset strip fins, perforated fins 
and louvered fins. A somewhat different concept for the reduction of thermal resistance 
and consequent enhancement in heat transfer is to induce longitudinal vortices in the 
flow field. Longitudinal vortices generated by delta wings or winglets attached to the 
surface at an angle of attack persist over a long streamwise distance. The spiraling flow 
of these vortices serves to entrain fluid from their outside into their core. Thus, these 
vortices disrupt the growth of the boundary layer and serve ultimately to bring about 
enhancement of heat transfer between the fluid and its neighboring surfaces at the price 
of relatively less increase in i:)ressui’e penalty. 

The geometrical configuration considered in this thesis are representative of single 
elements of either a compact gas-liquid fin-tube crossflow heat exchanger or a plate-fin 
heat exchanger. A novel numerical model has been developed to study the suitabil- 
ity of the proposed modification on the sm'face-geometry of the above mentioned heat 
exchangers. The model involves computation of three dimensional Navier-Stokes and 
energy equations in a rectangular channel with a built-in wing or winglet-pair. A de- 
tailed evaluation of the performance parameters with regard to enhancement of heat 
transfer and associated flow losses using various wing-type vortex generators has been 
accomplished. Prediction of heat transfer in the laminar regime is validated through 
comparison with experiment results of a research group in the Ruhr University, Bochum 
of Germany. The numerical model and the experiment corroborate each other reason- 
ably well. 

From heat transfer point of view, delta wing is found to be more effective than : 
winglet-pair. However, most convective heat transfer processes encounter two types 
of losses, viz., losses due to fluid friction and those due to heat transfer across finite 
temperature gradient. Because these two phenomena axe manifestations of irreversibil- 
ity, an evaluation of the augmentation techniques is also made from a thermodynamic 
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viewpoint. Conclusions that are drawn thus identify wiugk^t-pairs more efficient than 
the delta wings. 

The model is also extended to compute turbulent flows. The flow is described by the 
unsteady Reynolds averaged Navier-Stokes equations and the k-e model of turbulence. 
To mimic the transport processes near the wall, a two layer wall function treatment 
has been performed. The computed results are compared with the measurements of a 
research group in Stanford University of the USA. In general, the prediction is good 
although slight discrepancies are observed in the prediction of turbulent kinetic energy 
at places where the boundary layer is severely distorted by the vortices. Finally some 
important predictions are made with respect to enhancement of heat transfer for various 
Reynolds nmnbers in the turbulent regime. 

It is understood that the wall function approach used in turl)uk'.nt flow computations 
greatly reduces the storage and computational time as compared to those associated with 
the low Reynolds-number near- wall model of Jones and Launder (1973). However, it 
is also well-known that the wall-function approach is used on a rather simple intuition 
and it is not always consistent with complex physical picture near the wall. Instead 
of computing with the wall-function approach, two different but representative low- 
Reynolds number models designed for use with a standard k-e two-ecpiation turbulence 
model should be attempted. The models due to Lam and Brernhorst (1981) and Rodi 
and Mansour (1993) may be chosen for this pinpose. 
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Appendix A 

The Program Substructure 


In order to solve the full Navier-Stokes and energy equations for both laminar and 
turbulent flows, a computer code has been developed based on the numerical methods 
described in Chapter-3. The code has been written in FORTRAN-77 language. The 
code consists of several subroutines each of which has a set of specific tasks to carry 
out. The main program (DELTAM) and the subroutines are written in the form of 
modules. A flow chart is provided in Fig. A-1 showing the operational sequence of the 
various subroutines with their major communication links with the main program. A 
short description of different indices and the main functions of different subroutines are 
given below: 

DSriT 

IREST 


START 

RESTAR 


This is a subroutine which initializes the entire calculations for a set of 
input variables such as Reynolds number, Prandtl number, grid sizes etc. 

This is an index which denotes whether a restart from a set of existing 
field variables are desired (=1) or the computation should start from the 
beginning (—0) 

This subroutine creates the guessed field for the dependent variables 
and initiates computation 

This is a subroutine where the field variables are equated to a data field 
which waLS existing beforehand 
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CONTI 

JSTAM 

BCC 

BCCST 

BCO 

CEQCP 

EPSI 

BCNS 

BCNSST 

VELALT 

ITUR 


This is a subroutine which controls the pressure- velocity iteration loop 

This is a parameter which denotes whether there is a stamping on the 
walls (=1) or not (=0) in the computation for laminar flows 

This subroutine deals with the boundary conditions for continuity 
equation 

This subroutine deals with the boundary conditions for continuity 
equation in the presence of stampings on the walls 

In this subroutine some of the geometrical parameters for the obstacle 
are defined 

This is a subroutine where divergence of velocity vector (D) is evaluated 
and the pressure correction is done 

A small value prescribed as the upper-bound for the velocity divergence 
in each cell 

This subroutine deals with the boundary conditions for Navier-Stokes 
equations 

This subroutine deals with the boundary conditions for Navier-Stokes 
equations in the presence of stampings on the walls 

This is a subrotitiue which stores the current solution (velocity) for 
time level (n) in order to start a new time cycle of (u-hl) level 

This is an index which denotes whether a laminar flow(=0) or a turbulent 
flow (=1) is being computed 
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TURBU 

BCKD 

KDALT 

KDEQCP 

TURVIS 

DIVKD 

EPKD 

TICORR 

DIVMAX 

STAT 

ITA 


This subroutiiK' controls tlu' subroutiiK's for computation of tui-bulent 
viscosity witli tlu' luilp of k and e equations 

This subroutine deals witli the boundary conditions for k and c 


This subroutine stores the current variables of k and e in order to start 
next iteration 


This is a subroutiiK' where the k and e for the next iteration are 
explicitly ('valuat('d 


This is a sul[)routine where the turbulent kinematic viscosity for each 
cell is calculated 


This paraiucder denotes maximmn discrepancy of k and e between two 
iterative steps 

A small value that determines the convergence for the calculation of 
the equations for k and e 

In tins subroutine where the value of the time increment (At) and the 
dollar cell coefficient (op) are calculated 

This is a parameter which denotes the maximum time increment of 
field variables 

A small value which determines the conditions for steady solution 

Number of iterations in the time (one way) coordinate or number of 
time steps 


ITAMAX Maximum number of allowable iterations in the time direction 
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NSEQCP 

BCOV 

TIGRAD 

ENERGY 

BCT 

BCTST 

BOOT 

TEQCP 

DIVT 

TSTAT 


This is a subroutine where the velocities for the next time step are 
explicitly evaluated from the discretized Navier-Stokes equations 

In this subroutine the velocity boundary conditions for the obstacle are 
specified. This is used in conjunction with BCO 

This is a subroutine where the time-gradients of the velocity components 
are calculated 

This is a subroutine which controls the solution of energy equation by 
a SOR or equivalent technique 

This subroutine deals with the temperature boundary conditions for 
the confining surfaces of the channel 

This subroutine deals with the temperature boundary conditions for 
the confining surfaces in the presence of stampings on the walls 

This subroutine deals with temperature boundary conditions for the 
obstacle 

This is a subroutine where the temperatme field is calculated 

This parameter denotes inciximmu discrepancy of the temperature 
between two iterative steps 

A small value that determines the ct)nvergence for stc.'ady solution of 
the energy equation 


OUTPUT 


This is a subroutine which prints the output in the assigned output files 
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Fig. A-i 
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Fig. A- I 














